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NOMENCLATURE 


A Pre-exponential factor 
C Specific heat 
C Specific heat of unpyrolyzed wood 
C Specific heat of char 
d Quenching distance 
D diffusion coefficient 
2 2 
Ap Ch 
a a 


D Damkohler number, D = e 
r 


-£0 
E ~ Activation energy 


F Externally applied heat flux 


h Boundary layer thickness 


Convective heat transfer coefficient 


Sal 


Le Lewis number, Le = AP ED 
m Mass flux, m = pv 
M Nondimensional mass flux, M = meen 
q Heat release 
Q Nondimensional heat release, Q = SUAS - Te} 
R Ideal gas constant 
alse.) 
R Nondimensional reaction term, R = DY Y-exp[/——————] 
Ont ° 

l - a (1-6) 

tc Time 


uh Temperature 


v Velocity 
x Position coordinate 


ae Mass fraction 
Greek 


a Nondimensional Factor, Peed ace TL/Tp 

B Nondimensional activation energy, B = one 

B Nondimensional activation energy, B. = E/RT, 
6 Laminar flame thickness, 6 ~ Dr a se 
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6 Nondimensional temperature, @ = (T - T/T, - oe) 
aN Thermal conductivity 
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E Nondimensional position coordinate, € = x/h 
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Subscripts 
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CHAPTER 1 


INTRODUCTION 


The burning of wood and other cellulosic materials has long been 
of central importance to fire research because wood is a primary 
material for building construction. Therefore, a fundamental study to 
understand the chemical and physical processes that occur during the 
combustion of cellulosic materials is a significant step toward 


achieving effective control and prevention of unwanted fires. 


Cellulosic materials (such as wood, cotton, paper, etc.), when 
subjected to heating with sufficient intensity and for a long enough 
duration, will ignite to yield sustained combustion. Depending upon 
whether the ignition occurs with or without the aid of an external 
ignition source, the result is accordingly classified as spontaneous or 


piloted ignition. 


The objective of this thesis is to study the piloted ignition of 
wood theoretically. Although there are some theoretical models for the 
auto-ignition of pyrolyzing materials (further discussion in chapter 5), 
there is no theoretical model for piloted ignition of wood. Piloted 
ignition, like auto-ignition, includes the solid parHoleet g and gas-phase 
chemical reactions. However, the mechanisms for piloted ignition are 


different from those of auto-ignition. Actually, the detailed mechanisms 


Z 
are still unclear. From the experimental observations of A. Atreya 
(1983), the important factors for piloted ignition include: 
(1) pyrolysis of wood, which produces gaseous-fuels (such as CO, CH, 
and higher-molecular-weight hydrocarbons) and inerts (such as CO,, H,0). 
(2) diffusion and mixing of the fuel with the ambient air, which 
produces a flammable gas in the boundary layer. (3) pre-mixed flame 
ignition by a pilot ignition source. (4) quenching by the surface 
and (5) establishment of a sustained diffusion flame after the pre-mixed 


flames has vanished. 


This thesis develops a theoretical model which attempts to 
explain and to quantify the above experimental observations. The piloted 
ignition criterion used in the theoretical model is defined as the 
development of a sustained diffusion flame after the pre-mixed flash. 
This ignition criterion is in sharp contrast with ignition criteria used 
previously, such as: minimum char depth, minimum surface temperature and 
minimum fuel flow rate. It is important to emphasize that the ignition 
criterion used in this work is not an arbitrary threshold. Instead, it 
is a natural consequence of the equation solved in this work. This model 
would be helpful for understanding the important parameters of fire 
safety (such as increasing activation energy of wood, decreasing pre- 


exponental factor of wood and decreasing thermal conductivity of wood). 


In order to solve the mathematical model which includes both the 
pre-mixed and diffusion flames, a new numerical scheme has been 
developed. This method is herein called "a combined analytical and 


numerical solution method for chemically reacting flow"; the method is 


a 
described in detail in Appendix A. Usually, the numerical calculation 
of a pre-mixed flame requires very small time steps and grid spacing, 
but diffusion flame calculations can use a much larger time step, since 
the changes are slower. Consequently, the numerical calculations for 
the transient period spanning both pre-mixed flame and diffusion 
flame requires is very difficult, often causing over-flow for the 
numerical calculation. The combined analytical numerical method is 


comparatively stable during this transient period. 


In chapter two we present a brief literature review for the 
development of the numerical scheme for pre-mixed flame propagation, 


quenching, and pyrolysis and ignition of cellulosic solids. 


In chapter three a simplified piloted ignition model for the gas 
phase is developed. Here, irregular grid spacing is used for the 
numerical calculation. The finite difference expressions are the same as 
equation 3-100 of Anderson (1984). The computer program has been 
validated by the numerical calculation of the pre-mixed flame velocity, 
which is well known for CH,-air mixtures. Also, the program has 
been checked by calculating the minimum fuel flow rate for the steady 
state diffusion flame in two ways, and the diffusion flame location in 
three ways. In this model, the solid phase is assumed to provide a 
constant fuel mass flux at a fixed surface temperature. (In the real 
case, both heat flux and surface temperature vary with time.) Although 
this model is simplified, it reveals the basic structure of the piloted 
ignition process, which includes a pre-mixed flame in the earliest 


stage of piloted ignition; this flame may be quenched if the fuel flow 


A. 
rate is too low (lower than required to establish a steady-state 
diffusion flame) or if the ignition source is too close to the surface 
(which is the quenching of the pre-mixed flame). This model also shows 
that as the pre-mixed flames vanish a diffusion flame appears. This 
investigation also confirms that piloted ignition is related to the 
extinction of a diffusion flame. The simplified model of chapter three 


has been published in Combustion and Flame [Tzeng et.al. (1990)]. 


In chapter four the extinction of a steady diffusion flame is 
investigated both numerically and analytically. The analytical solution 
uses large activation energy asymptotics (AEA) to solve for the 


extinction limit of the steady-state diffusion flame. 


In chapter 5 a complete simplified piloted ignition model is 
developed which includes the details of the solid-phase pyrolysis. In 
this model, the effect of fuel concentration on the ignition delay time, 
ignition fuel flow rate, and ignition surface temperature is 
investigated. It is found that the assumption of 25% fuel and 75% 
inert (by mass) gives the most accurate prediction of ignition 
time, denaefon fuel flow rate and ignition surface temperature. A 
parameter study is also conducted, which shows that the activation 
energy of pyrolysis, pre-exponential factor and thermal conductivity of 
wood are very important factors for the theoretical model. Finally, 
chapter 6 summarizes the overall conclusions of this thesis and suggests 


various possible avenues of future research. 


CHAPTER 2 


LITERATURE REVIEW 


As already discussed, the model developed in this thesis must 
include descriptions of the thermal decomposition of the solid to 
produce fuel gases, mixing of fuel and air in the boundary layer, pre- 
mixed flame propagation originating from the ignition source, quenching 
of this pre-mixed flame by heat losses to the surface,. and establishment 
of a sustained diffusion flame in the boundary layer. [Io the author's 
knowledge, there are no piloted ignition models in the literature, Thus, 
this review basically includes discussions of (1) gas phase chemistry, 
and (2) solid phase pyrolysis. The piloted ignition model essentially 


is a combination of the above two phenomena. 


2.1 Gas Phase Chemical Reaction 


This part includes ignition, propagation, and quenching of a 
pre-mixed flame. The laminar pre-mixed flame speed and flame structure 
have been extensively studied in the past. Spalding [1956] first 
introduced a time-dependent numerical scheme for the solution of the 
equations for conservation of mass, energy and different chemical 
species. He assumed arbitrary initial profiles and then solved the 


equations using a standard finite-difference method. Difficulties 


6 
in the numerical procedure arise because the thickness of the 
pre-mixed flame is very small (about 0.1~0.5 mm, William 1985), thus a 
grid spacing much smaller than the flame thickness must be employed. 
Also, the chemical reaction rate in the governing equations is much 
larger than the diffusion rate, making the profiles very steep and 
forcing the use of small time cone Thus, numerical calculations 
of pre-mixed flame propagation requires both small grid spacing 


and small time steps. 


A numerical study of laminar flame quenching was performed by 
Aly and Hermance (1981); they found that the quenching Peclet number Pe, 
increases with decreasing fuel concentration. The theoretical results 


are in fair agreement with the available experimental data. 


Spalding (1953) found that when a combustible mixture is 
ignited by heating a slab of gas, the amount of energy added to the 
gas must exceed a minimum value for ignition to occur. Numerical 
integration for slabs of various thickness, initially raised to the 
adiabatic flame temperature, shows that ignition does not occur for 
slabs thinner than the pre-mixed flame thickness. For thicker slabs, 
however, a propagating laminar flame develops. This observation is 
important for numerical simulation of a pilot flame or an ignition 


source. 


2.2 Solid Phase Pyrolysis 


Combustible solids can generally be categorised as either 
(1) those which decompose volumetrically, or (2) those which decompose 
superficially. Wood belongs to the first group, while some polymers 
(such as PMMA) belong to the second group. In this study, only wood 


has been considered. 


Physically, the properties of wood are highly dependent upon its 
microscopic structure. Wood consists of cells or fibers, whose diameters 
vary from 0.02 to 0.5 millimeters. The length of the cell fibers ranges 
from 0.5 to 1.5 mm. The porosity (ratio of the volume of pores to the 
volume occupied by the cell walls) of real woods lies somewhere in the 


range 40-75%. (Kanury 1970). 


Chemically, wood is composed of three major constituients: (1) 
cellulose (50% by weight), (2) hemicellulose (25%) and lignin (25%) 
(Stamm 1964). Pyrolysis of cellulose occurs in two distinct path ways 
(Shafizadeh 1981). At low temperatures (200-280°C) dehydration of 
cellulose produces dehydrocellulose and water. This process is slightly 
endothermic (except in the presence of oxygen) and leads to the 
formation of char, water, and volatile gases such as CO,, CO, and 
hydrocarbons. The gases evolved in the dehydrocellulose route are 
primarily noncombustible and the char which remains can oxidize through 
surface combustion. At higher temperatures (280-340 C) endothermic 


depolymerization of cellulose leads to the formation of tar-like 


products which are highly condensible and constitute the main gaseous 


fuel to support a gas-phase flame. 


Pyrolysis of thick samples of wood involves both the physics of 
heat and mass transfer and the kinetics of chemical decomposition. 
Consider a thick slab of wood, initially at room temperature, exposed to 
a high-intensity thermal radiation source. The temperature of the solid 
gradually increases with time, the surface temperature being the 
highest. Before the surface layer decomposes, evaporation of moisture 
occurs and a moisture evaporation zone begins to travel into the solid. 
At later times, the pyrolysis zone begins to develop and then to 
propagate slowly into the interior of the solid, leaving behind a 


thermally insulating layer of char. 


Regarding the mathematical model of pyrolysis, the formation and 
growth of a char layer which protects the decomposition zone and the 
virgin material involves many physical and chemical processes, which 
makes it very difficult to develop a numerical solution that contains 
all of them. A few mathematical models of wood pyrolysis have been put 
forward, starting with the work of Bamford et.al.(1946). They treated 
wood as a solid of constant thermal diffusivity. Kung(1972) assumed the 
thermal properties of the solid varied continuously from their values 
for virgin wood to their values for char. Kansa et.al. (1977) included 
a momentum equation for the motion of gases relative to the solid and 
obtained good agreement with Lee et.al's (1976) experimental data for 
low heat fluxes. At high heat fluxes, however, poor atreement was 


obtained; this was attributed to ignoring the effects of structural 


changes (shrinkage and cracking) and to the assumption of a single 


step-pyrolysis reaction. 


An analytical solution for surface temperature was derived by 
Atreya (1983) using the integral method. The solution shows excellent 
agreement with the experimental results. Base on this surface | 
temperature, an analytical solution for the pyrolysis mass flux of wood 
has been derived by Atreya and Wichman (1989). This equation is obtained 
by assuming constant wood density. Thus, the solution is only valid in 
the early stages of pyrolysis. In this research, both the analytical 
(integral) solution and the numerical solution of the simplified wood 


pyrolysis equations have been used to describe the solid phase. 


CHAPTER 3 


A ONE-DIMENSIONAL MODEL OF PILOT IGNITION 


In this chapter a simplified model of piloted ignition is 
analyzed. The two-dimensional coupled solid and gas phase problem is 
simplified by assuming that the mass evolution rate from the combustible 
solid is a constant , and by employing a plane rather than a point 
ignition source. With these assumption the model problem requires only a 
transient one-dimensional analysis of the gas phase. The model equations 
are solved numerically using the fast scheme discussed in chapter 1 and 
Appendix A. The pilot flame is modeled as a thin slab of gas that is 
periodically heated to the adiabatic flame temperature of the 
stoichiometric mixture. The effects of : (i) the location of the 
ignition source, (ii) the fuel mass evolution rate from the surface, 
and (iii) the surface temperature of the solid are investigated. An 
explanation is produced for the pre-ignition flashes that are observed 
experimentally. A criterion for positioning of the pilot flame is 
proposed. The minimum fuel evolution rate, by itself, is found 
insufficient for predicting the onset of piloted ignition; heat losses 
to the surface plays an important role. Also, the conditions at 
extinction of a steady diffusion flame are found to be practically 


identical to those for piloted ignition. 
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3.1 Background 


The phenomenon of piloted ignition is rather poorly understood. 
This is evident from the numerous empirical ignition criteria that 
have been proposed in the literature. Some of these are: critical 
surface temperature, critical fuel mass flux, critical char depth and 
critical mean solid temperature. Of these, critical fuel mass flux at 
ignition appears physically the most reasonable, but critical surface 
temperature has proved to be the most useful, since it can be easily 


related to flame spread. 


The physical mechanism of piloted ignition is quite complex. The 
solid must first chemically Alert D to inject fuel gases into the 
surrounding air. This produces a flammable mixture (in the boundary 
layer), which is ignited by an ignition source. A plausible physical 
configuration for this process is illustrated in Figure 3-1.In the 
early stages of solid pyrolysis, the fuel flow rate is very small and 
the fuel-air mixture in the boundary layer is not combustible. As the 
fuel evolution rate increasses with time this mixture becomes rich 
enough to allow a premixed flame to propagate through the boundary 
layer. This premixed flame consumes nearly all the available fuel and 
is quenched by heat loss to the surface, unless the fuel evolution 
rate is large enough to replenish the consumed fuel. Thus, either a 
flash (a quenched premixed flame) or a sustained diffusion flame 
(piloted ignition) is observed. Experimental observations of this 
phenomenon are shown in Figure 3-2. The momentary rise in surface 


temperature prior to sustained ignition is due the flashes. 
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From experimental obserations, it is clear that the piloted 
ignition process involves both premixed and diffusion flames. Also, the 
true criterion for sustained piloted ignition is determined once the 
conditions that allow the conversion of a premixed flame to a diffusion 
flame are known. These conditions include magnitudes of the fuel 
evolution rate and the surface temperature that determines the heat 
losses responsible for quenching the premixed flame. Thus, a complete 
theoretical model of the piloted ignition process must include transient 
analysis of : (i) thermal decomposition of the solid to produce fuel 
gases, (ii) mixing of fuel and air in the two-dimensional boundary 
layer, (iii) premixed flame propagation originating from the ignition 
source, (iv) quenching of this premixed flame by heat losses to the 
surface, and (v) establishment of a sustained diffusion flame in the 
boundary layer. This rather formidable problem is considerably 
simplified by assuming that the solid-phase thermal response (to 
the applied heat flux) is a known function of time, and by employing 
a plane rather than a point ignition source. These assumptions reduce 
the problem to a transient one-dimensional analysis of the gas-phase 
phenomena, thus capturing the bare essentials of the piloted ignition 


process. 


3.2 Mathematical Model of One-Dimensional Piloted Ignition 


A diagram of the model problem is shown in Figure 3-3. A gas 
stream containing the fuel is fed uniformly and at a constant rate 


through a porous plate whose surface temperature is assumed constant. 
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Conditions at a distance h above the plate are maintained constant by a 
fast flowing oxidizer stream. A plane ignition source (whose distance 
from the plate, 0 < x. < h, is variable) is placed near the porous 
plate. This ignition source is periodcally "turned on" to test for 
piloted ignition. The chemical reaction is assumed to be a simple 


second-order, one-step irreversible Arrhenius type, 
Fuel + v 0, > Products=sdmtheat) (3-1) 


The governing energy and species equations are: 


OT MS Oost Amano. © ee 
p— + pv— = —|— —] + ae ve 87, (3-2) 
ét Ox @x|C_ ax C 
Pp P 
ax ay ) oY 2 
flees + pv— = — Lia's - Ap eet (3-3) 
at ox Ox Ox 
and 
oY oY a oY 2 
—— + pv— = — oD - vAp rohan (3-4) 
at Ox Ox Ox 
and the initial and boundary condition are: 
at to 0 toe ae 
Ye = 0, Yo = Ne ance is = he (3-5) 


L? 


oY, 
er ee) ee |)? pD— = pv(Y, - Probe 
Ox 
oye 
pO ea (256) 
° 
Ox 
T=T 
s 
Bioeeatexccmni. tO; 
Yp = 0; Yo = hare T = hae (3-7) 


To simplify these equations, p is assumed constant, implying that m = pv 
is not a function of x; thus, mass conservation is automatically 
Satisfied. Also, 4d, a; and D are assumed constants. 


The above equations are non-dimensionalized as follows: 


xX At m C h 
£ =|— ,-r = 54 ite 
h Coh » 
4 P 
2 2 
I< c. A pt .Gen if 
G@ = —— J, D= Lot Dees EXP[-B ] (3-8) 
T, - Ty r 
E q 
fe) ommne) 
Cee A ee he fe 
RT ae - cee) 
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and 


-B(1L - @) 
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the governing equations may be rewritten as 


6 Q 
i Y¢ = -1 ie (3-9) 
ye -vV 


The initial and boundary conditions become: 
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Equations (3-9) along with the initial and boundary conditions are 


solved numerically by using a finite difference formulation. 


3.3 Numerical Solution 


A nonuniform grid is used for numerical computations. As shown 


in Figure 3-4, this grid is'closely spaced around the ignition source to 
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ensure satisfactory resolution for the premixed flame propagation stage. 
An efficient numerical method that splits the solution procedure into a 
homogeneous reaction (explosion) calculation and a subsequent convective 
diffusive calculation was developed. This numerical procedure is 


described in Appendix A. 


3.3.1 Premixed Flame Velocity Calculations 


First, the premixed flame velocity was calculated in order to: (i) 
adjust the input parameters, (ii) verify the numerical code, (iii) 
determine an appropriate grid spacing and time step, and (iv) determine 
the input energy needed for the ignition source. For these 
calculations, the fuel is assumed to be methane and essentially the 
input parameters given by Coffee et.al. (1983) were used. However, the 
reaction order used in Coffee et.al. (1983) is 3 as opposed to 2 for the 
present problem. Thus, the pre-exponential factor, A, was slightly 
modified to obtain the same steady-state premixed flame velocity in 
a stoichiometric mixture. These input parameters are E = 29.1 Kcal/mole, 
Ty = 2230 K, T, = 298 K, vy = 4, q = 11355 cal/g fuel, 

Cr = 0.323 cal/g K, A= 1.5 x 10°. cal/em K sec, p = 3.74 x 10 a/cem 
aN 9 | 

Le = —— = 1] and pA = 3.56 x 10 /sec. 
a 

The premixed flame velocity calculations were performed by using 
the numerical code developed to solve the governing Equations (3-9), 
with the initial and boundary conditions modified as follows: (i) 
Initially a stoichiometric mixture of methane and air is assumed to 


exist above the solid surface and the fuel flow rate from the surface 
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is set equal to zero (i.e, m= 0). Thus, Y¥_(0,x) = 0.0546 and 
Y¥5 (0.x) =0.219. (ii) The surface temperature is assumed equal to the 


adiabatic flame temperature, i.e. @ = a = 1. 


Thus, the premixed mixture is ignited at the surface and the flame 
propagates toward the outer boundary, eventually attaining a steady 
speed. This steady-state flame velocity is determined by measuring the 
propagation rate of the temperature profile and by assuming the 
densities of the burned and the unburned gas mixtures to be identical. 

A uniform grid was employed for these calculations and the influences 

of diffenent grid spacings and time steps were investigated. The results 
of these calculations are summarize in Table 3-1, which shows that the 
calculated steady-state flee velocity is very close to other 
calculations, as well as to the experimental result, thus giving 


credibility the numerical code and the input data. 


The grid spacing and the time step for the piloted ignition 
calculations is determined from Table 3-1, which shows that a grid 
Spacing of A€é = 3./ x 10 i provides sufficient resolution. This 
corresponds to an actual physical spacing of 0.055 mm, which is smaller 
than the flame thickness; this value is estimated by using 6 = EE 
Williams (1985), where Po and v, are the unburned mixture density and 
velocity respectively, giving a flame thickness of 0.092 mm. A finer 
grid was not used in order to reduce the CPU time. Since the premixed 
flame propagation speed is not the desired final result of this piloted 
ignition study, the largest time step shown in Table 3-1 RGR 
is used to save CPU time. Also, twice this time step (47 = 5.4 x 10°) 


was used prior to turning on the ignition source. 


oe 


3.3.2 Numerical Simulation of the Ignition Source 


In the numerical code, the ignition source is turned on by simply 
raising the temperature of a slab of gas to the adiabatic flame 
temperature. The thickness of this slab is chosen to be three grid 
spaces, which is approximately twice the premixed flame thickness. 
Choice of this thickness is based on the observation Williams (1985): 
"Ignition will occur only if enough energy is added to the gas to heat 
a slab about as thick as a steadily propagating adiabatic laminar flame 
to the adiabatic flame temperature." The reason for choosing an 
ignition source thickness larger than the premixed flame thickness is 
that the gas mixture is not stoichiometric. Thus, the chosen thickness 
ensures ignition if the fuel/air mixture in the boundary layer has 


reached its lean flammability limit. 


It is important that the amount of energy added to the system by 
the ignition source be minimized. Thus, the pilot cannot be left on 
continuously. This is especially critical because a plane rather than 
a point ignition source is used to render the model problem one- 
dimensional. To overcome this difficulty a novel method has been 
devised. When the ignition source is turned on the data for temperature 
and species is placed into a temporary file. After calculating for 
several time steps the occurrence of ignition is determined by observing 
the development of the temperature and species profiles. If ignition is 
not observed, then the data from the temporary file is recalled and the 
calculation is restarted at that point. This procedure eliminates any 
accumulation of energy from periodic trials prior to the occurrence of 


piloted ignition. 
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3.3.3 Steady-State Diffusion Flame Calculations 


In the latter stages of piloted ignition, the premixed flame 
develops into a diffusion flame, which slowly moves toward its steady- 
state location. The governing equations and the boundary conditions 
describing this steady diffusion flame are simply the steady-state 
versions of Equations (3-9), (3-11) and (3-12). Thus, the location of 
the steady diffusion flame and the corresponding temperature and species 
concentrations may be found either from the large time solution of the 
time dependent Equations (soe (3-10) and, (3-12), (i.e... ini the later 
stages of piloted ignition) or by directly solving the steady-state 
versions of Equation (3-9),(3-11) and (3-12). It may also be obtained 
analytically in the thin flame sheet Burke-Schumann limit, which 


gives Kanury (1975). 


; geenyelevyas Yea, \ soe (213) 


Table 3-2 shows a comparison of the diffusion flame location 
calculated by all three methods. The transient piloted ignition 
calculations were terminated after about 4 seconds (elapsed real time) 
because changes in the diffusion flame location had become negligibly 
small. In these calculations the surface temperature of the solid was 
assumed to be 298 K and the ignition source was located at the minimum 
distance from the surface. (This minimum distance for the ignition 


source will be discussed in Sec. 3.4.4) 


Numerical solution of the steady-state versions of the Equations 


(3-9),(3-11) and (3-12) was obtained by iteration. Calculations were 
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started by guessing the location of the diffusion flame and the 
corresponding temperature and specied profiles. Using this initial 
guess, the non-linear reaction rate term (R; See Equation (3-9)) was 
calculated and then the energy and species equations were solved to 
obtain a new set of conditions. This iteration was performed until the 


solution converged to within a prescribed error. 


It is evident from Table 3-2 that the results of all three 
calculations compare quite satisfactorily. These calculations further 


confirm the validity of the numerical code for the piloted ignition 


calculations. 


3.4 Results and. Discussion 


After verifying the piloted ignition model for the limiting cases 
of premixed flame propagation (the early stage) and development of a 
steady-state diffusion flame (the later stage), several aspects of the 
piloted ignition process were investigated. In this investigation the 
fuel concentration in the solid phase was taken as unity (Ye, =1), 
and unless otherwise specified, the solid surface temperature was 


ts (T. = Te. = 298 K). Various case studies are summarized below. 


fhe 


3.4.1 Piloted Ignition Phenomena 


To investigate the details of the piloted ignition process, the 
_§ 2 
following conditions were assumed: h = 1.5 cm, m= 8.34 x 10 g/cm sec 


T. =ez90 K,°Y. = 1, and the ignition source is placed at 1.186 cm 


fs 
(which corresponds to the theoretically predicted location of the steady 
diffusion flame). Note that in experiments both surface temperature and 
fuel mass flux will increase with time as the solid is heated by a 
contant heat flux. Also, the fuel concentration in the solid is 
typically less than.100% <iver, Yeu < 1) because of the presence of 
moisture. Thus, the calculated time for piloted ignition cannot be 
compared with experimental values. However, the gas-phase ignition 
process (which is the focus of this study) is not expected to be 
different. Furthermore, in the numerical model, it is fairly 


straightforward to assign the experimentally measured values of m,T. and 


Ye. once they become available. 


Figures 3-5 and 3-6 show the results of piloted ignition 
calculations. The mixture was ignited at t = 1.089 seconds, as seen by 
the sharp temperature peak and corresponding dips in the fuel and 
oxidizer concentrations. The premixed flame then travels quickly into 
the unburnt mixture in both directions i.e., both toward and away from 
the porous plate. This can be seen more clearly from the heat flux 
profiles (which are proportional to the temperature gradient) show in 
Figure 3-6. At t = 1.102 seconds (i.e. 13 milliseconds later) the 
premixed flame has consumed nearly all the available oxygen on the fuel 


side and all the available fuel on the oxidizer side of the ignition 
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FIGURE 3-5 Histary of the piloted ignition process from ignition (f = 1.089 s) to development af 


a stcady diffusion flame (f = 3.887 s). Note: the scales for Ye and Y- arc nat shown for 
clarity. They differ from the scale for @ by a constant factor that is 0.250 far Yo and 0.167 
for Ye. 
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source, thereby establishing conditions appropriate for the formation 
of a diffusion flame. At this time the surface experiences a large 
heat flux (see curve (3) in Figure 3-6), which is probably responsible 
for the sharp momentary rise in temperature observed during the 
experiment shown in Figure 3-2. Note that during the premixed flame 
stage the heat flux is very large. However, once conversion to a 


diffusion flame begins the magnitude of the heat flux drops sharply. 


At times larger than 1.102 one the temperature, fuel and 
oxygen concentrations slowly adjust to those for a steady state 
diffusion flame, which is finally established at t = 4 seconds. Note 
that its final location is nearly identical to that of the ignition 


source (&£ = 0.28). 


3.4.2 Quenching of the Premixed Flame (Flashes 


If the ignition source is too close to the cold porous wall (or 
the sample surface) then thermal quenching of the premixed flame 
prevents the development of a diffusion flame and hence the occurrence 
of piloted ignition, resulting in a flash. Similar behavior is observed 
if the fuel flow rate is lower than the lean flammability limit. This 
is discussed in the next section. In the case discussed here, the fuel 
flow rate is the same as that for the case in Section 4.1, but the 
ignition source is placed closer to the surface (i.e., at x, = 0.193 em, 
as compared with x, = 1.186 cm for the previous case). For this case, 
gas ignites at 0.1 seconds (much smaller than 1.089 seconds for the 


previous case) but the ignition is not sustained. 
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Figure 3-7, 3-8, 3-9 and 3-10 show respectively the temperature, 
fuel concentration, oxygen concentration and the heat flux profiles 
during this momentary flaming. It can be seen that nearly all of the 
oxygen on the fuel side and the fuel on the oxygen side are consumed, 
but the temperature eventually decays and the fuel and oxygen are 


subsequently replenished. 


It is interesting to note that the fuel concentrations at the 
location of the ignition source for both this and the previous case are 
nearly the same (Y, = 0.0418 for the present case and Y¢ = 0.0414 for 
the previous case) and yet sustained ignition was not achieved here. 
Table 3-3 shows the fuel concentration at the location of the ignition 
source and at the time of piloted ignition for different fuel flow rates 
and ignition source locations. Note that the value of Y¢ is remarkablely 
constant (average is 0.0445) despite the fact that in many cases the 
flame was quenched. This value of Yp is approximately 80% of the 
stoichiometric mass fraction and is about 1.74 times the stated lean 


flammability limit for CH, [Kanury(1975)]. The reason for the difference 


4 
between the lean flammability limit and Ye at ignition is not clear 

to the authors. It may occur because of the assumption of one-step 
chemistry or because of the difference between the initial and boundary 
conditions. Further theoretical and numerical analysis of this 
phenomenon are obviously necessary. However, regardless of their outcome, 
it is clear i) that piloted ignition cannot be predicted solely on the 


basis of the lean flammability limit, and ii) that heat losses to the 


solid surface play a very important role. 
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3.4.3 Minimum Fuel Rate for Piloted Ignition 


In this case the fuel flow rate was reduced to 3.128 x ica ep 
sec, which is slightly less than the minimum fuel flow rate listed in 
Table 3-3 (3.34 x a acne sec) for which piloted ignition was 
possible. Note that this fuel flow rate is roughly one-third the fuel 
flow rate used in the previous two cases. Here, regardless of the 
location of the ignition source, a sustained diffusion flame was not 


obtained. 


Figure 3-11 through 3-14 show the temperature, fuel concentration, 
oxygen concentration and heat flux profiles during flashes caused by the 
low fuel flow rate. For the results presented in these figures, the 
ignition source was located at 0.67 cm from the surface, which 
corresponds to the theoretical location of the steady diffusion flame. 
It can be seen that although most of the oxygen on the fuel side and the 
fuel on the oxygen side are consumed by the premixed flame, the fuel 
supply rate is too low to enable the establishment of a diffusion flame. 
Hence, the temperature quickly falls to its ambient value. A comparison 
of Figure 3-8 and 3-12 shows how slowly the fuel is replenished. In 
Figure 3-8, the fuel concentration at the surface recovers to about 
75% of its initial value in 0.1 seconds, whereas in Figure 12 the fuel 
concentration at the surface recovers to only about 50% of its original 


value in 0.5 seconds. 


Since for successful piloted ignition it is necessary to establish 
a steady diffusion flame, the minimum fuel flow rate at extinction of 
the steady diffusion flame is expected to be close to that required 


for piloted ignition. The minimum fuel flow rate at extinction is 
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determined by numerically solving the steady-state version of the 
governing Equations (3-9). For this calculation the temperature and 
species profiles are first calculated for conditions under which the 
steady diffusion flame is known to exist. These results are then used as 
an initial guess for subsequent calculations with slightly lower fuel 
flow rate. This procedure is continued until the temperature and species 
profiles are the same as those with no chemical reaction, i.e., no 
diffusion flame exists. This rather stringent condition for extinction 


produces a lower bound for the fuel flow rate. 


By using the above procedure the minimum fuel flow rate at 
@x-iiec.0on 1s found to be 2 7667x 10 af Mpa sec. This fuel rate 
is only 8% lower than the minimum fuel flow rate for piloted ignition 
Cee GK 107° PiGy sec). This result seems to substantiate the 
hypothesis that conditions at extinction of a steady diffusion flame 
are similar to those at piloted ignition. This hypothesis was used by 


Atreya and Wichman (1987) to obtain an approximate analytical solution 


for the piloted ignition problem (solid-phase solution only). 


All results presented thus far assume that the surface temperature 
of the solid is held constant at ay irene. ay = a = 298 K). As T. eS 
increased, heat losses to the surface are decreased; this is expected to 
reduce the minimum fuel flow rate required for the establishment of a 
steady diffusion flame and for sustained piloted ignition. Calculations 
with T. = 491 K show that the minimum fuel flow rate at extinction of 
a diffusion flame reduces to 2.43 x Oe eee sec and that for 

5 2 


sustained piloted ignition reduces to 2.71 x 10 g/cm sec. Once again, 


the fuel flow rate at extinction of a diffusion flame is about 10% lower 
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than for sustained piloted ignition. 


3.4.4 LOCatCLON Of. Che I enitlon source 


As noted in Section 3.4.2 the location of the ignition source is 
very important in determining whether or not sustained piloted ignition 
will occur. This fact was experimentally observed by Simms (1963), who 
states: "The pilot ignition time depends not only upon the intensity of 
radiation and the density of wood, but also upon the position and 


possibly upon the size of the pilot flame as well." 


Table 3-3 summarizes the calculations for different locations of 
the ignition source at a given fuel flow rate. It can be seen that there 
exists a minimum location of the ignition source for sustained piloted 
ignition regardless of the fuel flow rate. For ignition source distances 
lower than the minimum location, the flame is quenched (see Section 
3.4.2). As the fuel flow rate is decreased, the minimum location 
approaches the theoretical location of the steady-state diffusion flame. 
Thus, the optimum location of the pilot is the location of the steady- 
state diffusion flame, which can be estimated from Equation (3-13). 

Note that in the experiments, the fuel flow rate and surface temperature 
slowly increases as the solid is exposed to a given external radiation. 
Thus, to find the minimum time (or fuel flow rate) at which sustained 


ignition occurs, the pilot must be located at the eventual position of 


steady diffusion flame. 


As the surface temperature increases both the minimum fuel flow 


rate for sustained piloted ignition and the minimum distance of the 
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ignition source decreases. A plot of the minimum distance of the 
ignition source versus fuel flow rate at surface temperatures of 298 K 
and 491 K is shown in Figure 3-15. For T. = 491 K the minimum distance 
is generally much lower than for T. = 298 K. However, at the higher 
fuel flow rates both curves approach the Same vertical asymptote, 
x = 0.18cm. It is interesting to compare this value to the 
theoretically predicted quenching distance for premixed flame 
propagation. Thus (see Williams (1985). pp. 268-276, and Sec. 3.3.1), 
co = ab = OILERS, = 40(0.092 mm) = 0.37cm, which is twice the value 
of x; found here. This is reasonable, because the ignition source is 
periodically raised to the adiabatic flame temperature while quenching 
temperatures are usually much lower. In addition, the constant factor 
a may vary significantly with chioce of parameters (up to 50%,see 
Williams (1985); also,the derivation for So in Williams (1985) assumes a 
circular tube, not a plane wall. The horizontal asymptote in Figure 3-15 
suggests a minimum fuel flow rate, below which ignition of the gas is 
impossible; this may be interpreted as a blowoff limit due to thermal 


quenching. 


a. de oLeniticance 


The simple one-dimesional numerical model presented here reveals 
the basic structure of the piloted ignition process. In the early stages 
of piloted ignition a premixed flame quickly propagates through the 
unburned mixture, consuming nearly all the fuel on the oxidizer side and 
all the oxygen on the fuel side, thus establishing conditions that are 


appropriate for the formation of a diffusion flame. This premixed flame 
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may be quenched (resulting in a flash) if the fuel flow rate is too low 
or if the ignition source is too close to the surface (due to excessive 
heat losses). The analysis confirms that the sharp rises in surface 
temperature during experiemnts (flashes in Figure 3-2) occur because of 
burnout of the accumulated fuel during the premixed flame propagation 
stage. Note that during this stage the heat flux is an order of 


magnitude larger than in the steady diffusion flame stage. 


Under appropriate conditions a diffusion flame is established, 
which slowly moves toward the steady-state diffusion flame location. It 
has been found that the minimum fuel flow rate in itself is not 
sufficient to predict the onset of sustained piloted ignition (as is 
often cited in the literature) and that heat losses to the solid surface 
play an important role. These calculations also confirm the hypothesis 
that the conditions at extinction of a steady diffusion flame are very 


nearly identical to those required for piloted ignition. 


The optimum location of the pilot flame for actual experiments is 
found to be the eventual location of the steady diffusion flame. As 
expected, both the optimum location and the minimum fuel flow rate 
required for piloted ignition decrease with increase in surface 
temperature. In experiments where the surface temperature increases with 
the time of exposure to external radiation, an incorrect placement of 
the ignition source will eventually result in ignition, but at a higher 
surface temperature and after a larger ignition delay time. This may 
partially explain the wide range (300 to 540°C) of measured surface 


temperatures for piloted ignition of wood. 
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TABLE 3-1: Flame velocities calculated by using different grid 
spacing and time steps. (units: cm/sec) 


NG 
5 ft0%* es aeTo.? 
ME 


NOTE: Coffee’s (1983) calculations--39.8cm/sec. 
Estimated from analytical formula of Williams (1985)--39.18 cm/sec. 
Experimental measurements from Kanury (1975)--37.3 cm/sec. 


45 


TABLE 3-2: Steady state diffusion flame location (x,) 


Fuel flow Xe from x, from x, from Flame temperature 
rate transient| steady- Equation | Numerical | Analytical 


2 calcula- state cal- (339 ) CPK) ( K) 
(g/cm sec)]| tions culation analytical 
after 4 (cm) (cm) 
sec(cm) 


225x10° 


8.34x10> 


4.17500" 


375x110" 


3.34x10> 
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TABLE Go =o. Location of the ignition source and the fuel 
concentration at this location. All calculations 
are Lor da = ih = 298 K 


Fuel flow rate| Location of | Y, at the |Time elapaed} Observations 
2 the ignition] ignition |for piloted 
(g/eom sec) source (cm) source ignition 
location . (sec) 
and at the 
time of 
ignition 


ignition 
ignition 
ignition 
ignition 
flame quenched 


oOo Oo © OF 
FREE 


ignition 
ignition 
ignition 
flame quenched 
flame quenched 


0 4. 
0 4. 
0 4. 
0 4. 
0 4. 


ignition 
ignition 
ignition 
flame quenched 
flame quenched 


SOLO voreS 
FFE 


ignition 
ignition 
flame quenched 
flame quenched 
flame quenched 


See S&S 
Js iss des dss dS 


No sustained diffusion flame regardless of the 
location of the ignition source. 


+ Location of the steady-state diffusion flame. 
* Minimum location of the ignition source for sustained flaming (piloted 
Lens eon) 


CHAPTER 4 


EXTINCTION OF A STEADY STATE DIFFUSION FLAME 


In the later stages of piloted ignition, the pre-mixed flames 
vanish, leaving behind a diffusion flame. This diffusion flame then 
moves slowly toward its steady-state location. As the flame reaches 
the steady-state condition, i.e., when all derivatives with respect to 


time disappear,the governing equations (3-9) become: 


Energy: 
2 
d @ d 6 
Me ee a OR, (4-1) 
dg 4d€ 
Fuel: 
a Y 
dy 
eet = a and se (4-2) 
Gee di § 
Oxygen: 
2 
dy d Yo 
M—- = ——> - ve. (4-3) 
dé d € 


In this chapter, the governing equations for the steady-state diffusion 
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flame have been solved numerically and analytically. 


The structure of one-dimensional flames was successfully 
analyzed by Linan (1974) by using large activation energy asymptotics 
(AEA) for the model problem of a counterflow diffusion flame. The 
boundary conditions for the steady state diffusion flame in this thesis 


are different from Linan'’s model, whereas the governing equations are 


very similar. 


In section 4-1, the equations are solved by using the finite 
difference numerical scheme, while in section 4-2, the governing 


equations are solved using AEA. 


4.1 Numerical Solution for the Steady-State Diffusion Flame 


The steady-state governing equations of the diffusion flame can 
be solved numerically. A finite difference method with iteration of 


the non-linear term was used. The finite difference representations are 


as follows: 


de /dgé=(6, 4-9; )/dE, 


2 2 2 
dng /de =(G,,°20,+6,_5)/d§ : 
Similar definition are applied to the fuel and oxygen concentrations. 


The finite difference equations are non-linear and conjugated. An 
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iteration method has been used to solve these equations. This procedure 


is described below: 


(i) Assume that there exists a diffusion flame, and that the 
temperature and species profiles are smooth near the 
reaction zone. 

(ii) Use the initial conditions in (i) to calculate the non-linear 
reaction rate term R, and then solve the temperature 
equations using the Thomas algorithm. 

(iii) Assume that Y- is the only unknown. Then the equations for 
the fuel concentrations are linear; again, by using the Thomas 
algorithm, the equations are easily solved. 

(iv) Now assume that Yo is the only unknown, and solve the oxygen 

equations. 
(v) Compare the solutions to the initial conditions, if the 
_4 
maximum error larger than a prescribed value (say, 10 ), 
then re-define the new initial conditions as: 
(New initial conditions) = 0.7*(old initial condition) + 0.3*( 
the, solutionsacalculated insteps.(ii),,,Ciii) «andy (ivo)s 
(vi) Use the new initial conditions to calculate the non-linear 


term, then solve the temperature equations again and repeat the 


procedure. 


By using this numerical code, the minimum fuel flow rate for diffusion 


flame is also calculated according to the following algorithm: 
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(i) First, assume that the fuel flow rate is high, and calculate 
the temperature and species profiles by using the steady- 
state numerical code. 

(ii) Using the solution in (i) as a new initial condition, then 
reduce the fuel flow rate slowly, and recalculate the 
temperature and species profiles. 

(iii) Continue step (ii) until there is no steady-state diffusion 
flame.(ie., the temperature and species profiles are the same 


as without chemical reaction). 


4.2 The Extinction Damkohler Number for Diffusion Flame. 


The (AEA) method used here has been successfully applied to 
other diffusion flame problems as well. C. K. Law (1975) solved the 
ignition and extinction of droplet burning by following a procedure 
similar to Linan’s (1974). Puri and Seshadri (1986) estimated the 
activation energy and pre-exponential factor for counterflow diffusion 
flames by using the extinction Damkohler number of a diffusion flame. 
The physical configuration for piloted ignition is different from 
these cases. The objective here is to calculate the minimum fuel flow 
rate for sustaining the one-dimensional flame sheet by using (AEA) 


(similar to Linan (1974)). 


It is well known that the ignition and extinction Damkohler 
numbers for a diffusion flame can be described by an S-shaped curve 


[Williams (1985)]. In such a curve, for Damkohler numbers smaller than 
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the extinction Damkohler number, having a diffusion flame is impossible 
(see Figure 4-2). Likewise in the numerical simulation of piloted 
ignition, when the fuel flow rate is too small, no sustained piloted 
ignition is possible. Thus, it seems reasonable to assume that sustained 
piloted ignition can occur only if the Damkohler number for the 
corresponding boundary layer is at least as large as the extinction 


Damkohler number, which is determined in the following analysis. 


The species and energy equations for the steady state diffusion 


flame are rewritten as follows: 


ye T q/C 
d d Pa aed Sas 
a ee en eee lek = | 1 |AAaY NX EXPY-E/RT] (4-4) 
ax dx dx 
pD Yo gs) ae 


The boundary conditions are: 


At x=0: 
dy. 
= pe ae m(Y 5. 5 Ye), 
dx 
dy, 
plo" = mY | (for no diffusion flame), 
dx 
Yo = 0 (for a diffusion flame), 
TL = T,(t), 
where T. is the temperature of the lower surface. 
At x=h: 


=P". 


Dicer! lees a Lime 
u 


co 


where Ty is the temperature for upper surface. 


These governing equations are non-dimensionalized with the 


A = te o 
following quantities; L = — (length), v = — (velocity), Lae 
— p a 


(temperature), where Yp a is dependent on the fuel flow rate. 


i. d 
Assuming pv = m= constant and Le = —=1 (— = pd), 
C_pD C 
P Pp 
and also defining ese PPG 8 ae ATS” x = x/L = x F : 
pr A 
and a= Fo SEC) UE VpAY . as convection time/reaction time, 
Cm : 
P 
¥ E 
et 7 land Tv. =-— gives 
RT, 
d d e if 
(Gath pee? |A e a | acak el ey eae (4-5) 
d x dx 
Yo ail 
The boundary conditions are: 
At x=0: be = 0 (for a diffusion flame), 
oe any ie (for no diffusion flame), 
iP pe A, 
T 
pat --E, 
- 


Fe) 


At x= 2 = Meee vrs (Oe 


(where Nhs eae We Ra RS 


Here, Yo a is dependent on the fuel flow rate with no diffusion 


flame. To eliminate the convection term (in order to eliminate the 


x 
e 


convection term, \ must be a constant), define ¢ = 1 - >: Less 
e 


obvious that at x = Os0C va ie eee Sand at x = 2, ¢ = 0. Thus, we 


have 


age pe I 
d Dy -y EXP[-T_/T] 
—z|y,]- | 1 |—=2*—+. (4-6) 
= (1/p,-(1-2)) 
Ye 1 


where z = (Ores Ds Sina The boundary conditions are now 


At x=0, (z=0): y,70 (for a diffusion flame), 
ya) an (for no diffusion flame (dependent on m)), 
ye-l, St ere to be determined by m), 
T-T, . 


Since me > Di; let B = te - rT a Then, using a Schvab-Zeldovich 


formulation with 


S (4-7) 
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re 2 - =: 
gives d Bb, /dz - 0, i= 1, 2. The solutions foryp 7 and) fo are 


By = Tty, = ©,,+C212 = aztT, for a diffusion flame 


= aztT ty. 9 (t-2) for no diffusion flame (4-8) 


B. = T + bigfoe Cy, + Gooz = His Fe A 


The frozen flow temperature an (defined as the temperature when the 


chemical reaction vanishes) is given by 


T, = qT, + (TL . TDG - Zz) = T +9 BCL az) (4-9) 
Therefore, 
pon az + Ws aM for a diffusion flame 
(4-10) 
= az + Ty eT vn oft - z) for no diffusion flame 
and 
Pe ON uy ibe 2. (4-11) 


Three local flow regimes exist in the limit T. + o (infinite 


activation energy) 


(a) Frozen flow: 


——se= 0 (no reaction), 


(b) Equilibrium flow with voy Os 
T=T, + az, 


cm 


(c) Equilibrium flow with DS Os 


Fe 


For all regimes 
(a) Nearly frozen (Ignition) regime: 


Chemical reaction almost equal to zero. 


T= T, = T + BCL 2) 


= az + Ye 1 keoez) 


MS 


oh (4-12) 


wie Yo of 4 Yo,0 
Yer =. 1-z 


(b) Partial Burning: 


Two frozen flow regimes separated by a thin reaction zone. 


(see Figure hole 
(c) Pre-mixed Flame: 


Frozen flow and equilibrium regime separated by the reaction 


zone. 
(d) Diffusion Flame: 


Two equilibrium regimes separated by a thin reaction zone. 


(see Figure 4-2). 


4.2.1 Analysis: 


Near Equilibrium, Diffusion Flame Regime: 


Let Zo be the location of the diffusion flame. 


56 


} . 
wcananacacccannacaneusange ete ehaad CEagAME GE  CORGALL En tne tt eRe Re On ZC € 


Frozen flow 


Partial burning 


Tm 


ee 


— Frozen ; 
ol 


FIGURE 4-1 Partial burning. 
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FIGURE 4-2 Diffusion flame. 
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pfeve Anes zi, we have yaa Oe T= 7a + (a - B)z, wheras for z > Zo: 


we have 


Ae Of Ta= Ty. 91 Glet f)2 ae Ae 2a is 0 and 


TL + (a - B)z = TL +1-(1+ p)z|, a ze Therefore, 


1 Cn * Pea eae e 
ZS —, giving T. = qT + = Forez.- Ze OCEej- 
l+ea 


ee = 2 
(the reaction zone) we put T = T, - ek € - €k,B, - € Ko Bouse 


eb} 
LA Lhe 
and & = ———  B, where « = — << l, We first match toetue 


€ Ly 
a 


right of Zo: with the same slope for z > Z.: Let 


1 1+ 8B 
ern - k) =; 
k, B 
dp, 
os >) | em Leis asboundaryscondition: 
OS tn Soo 


Now matching the slope to the left of Zo» we let 


1 a- B 
en - k = 1; 
ky nm 
dp, 
: —— aes = -l is the other boundary condition. 
dé + -« 


Again matching the slopes to the right and to the left of Zo from 


Equation (4-6) we obtain: 


Bo 


ze=l 
re 
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Sa er) wet oe T= te 
yeas 
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FIGURE 4-3 Near equilibrium, diffusion flame. 
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2 eS ed 
a Bp, DS € 4k, exp[-T_/T,] es 4 ig 

=~ = 2 = £58, -) (8, 48 exp [- (k,€4k, By). 
dé [1-6 (@/(14a))] (14) 


1 
3 
Defining k,...=)5 (reduced Damkohler number) 


Oe ey or 
oO 


and 
noe * EXP | 17.) 
Se x 
5 - —____=;-—; (4-13) 
[1-¢ (@/(14+e))] (142) 

gives 

2 

dp, - 7 - ue 

an, wenn Sins) (Pi +5 )ERE[-(k -o4+k, 0, )] 

dé 6 +. €6, 

= (B,-€) (By +8) EXP[-6 7 °/ 7 (k 5° 7&48,) 1. 

Therefore, 

2 

are = ~ ivama ian’ 

—z— = (B, - €)(B, + €)EXP[-60 °° (6, +r€)], (oa 1) 

dg 

173 sigh : 

where r= ko6 6 : The boundary conditions are: 

dB, ~ 

>a u as E ig ES 

dé 

dB ze 

a -l1 as € fi Sets 

dg 


Linan’s (1974) analysis shows that when bo is too small, there 


is no solution for the above differential equation; the extinction 
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Damkohler number en is approximately given as 
2 3 4 
Gate = e((l-r)-(l-r) +0.26(1-r) +0.055(1-r) }. (4-15) 


In other words, when Das < 8 oe? there is no diffusion flame (or the 

flame is extinguished). If the activation energy and the pre-exponential 
factor are fixed (or a particular fuel is chosen for calculations), the 
reduced Damkohler number 6, and the extinction Damkohler number es 


are then functions of the boundary conditions. 


In the one-dimensional piloted ignition model which is of concern 
here, when the fuel flow rate too small, there is no sustained 
diffusion flame after the flash. This may occur because the reduced 


Damkohler number is too small. The calculations of the Damkohler number 


are presented below. 


4.2.2 Calculation of the Extinction Damkohler Number. 


The fuel used here is methane and the physical properties are the 


same as given by Coffee (1983). We put: 
9 
ph = 3256 exe lO Cl sec), oe = 0.323 cal/g K, h = 1.5cn, 
_5 ee 4 
A = 6:2. x°10)\ Ki /m Kisec “GQ00M Reiat mais epee nD On x 10. g cal/cm 
sec K, Aeon 11355 cal/g of fuel, E = 29.1 K cal/mole, R = 0.001983 


Kcal/mole K, E/R = 14673 K, and v = 4. For every M (M = ene) we have 


£.0 ee oH “és + bpm Adie where Yeo = 1 for pure fuel, and 
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cya 

r o~ Lak, o = 

pas ee an = 35155 Bays (K), T. = See T/T, 
£f,0 Pp * 


(DEE 29B5 K)., ner eae Dee inate athe = 2 


vent /t, 


a a- B pr us T 
T. + , aes a1 wi mv pAY 
- * lt+a Cm 


| 
I 


~ 
II 


1 - aeey Also 6. and ee are given by Equations (4-13) and (4-15) 


and we use T. = 298 K and Y =]. 
LL ts 


For a given fuel flow rate m, the reduced Damkohler number and the 
extinction Damkohler are then easily obtained. The extinction Damkohler 
number is calculated by employing the following procedure. 

(i) Assume a higher fuel flow rate, 
(ii) Calculate the corresponding values of 6. and b oe? 
(iii) Check if bo < ae if not, reduce the fuel flow rate 


and repeat step (ii). 


This procedure shows that when the fuel flow rate is smaller than 


vs 2 
Sey x 10 p/om sec, 6 isesmaller*than 6  , (6 =0.827, 6 =0.86). 
re) oe fe) oe 


4.3. Comparison of Results 


Here the results of minimum fuel flow rate for piloted ignition are 


compared with the steady state numerical model and the analytical 
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Damkohler number analysis. The following results are obtained. 


shy, 


When the fuel flow rate is smaller than 3.13 x ie sient sec, 
there is no sustained diffusion flame after the flashes. This 
result is obtained by numerically solving the piloted ignition 
equation (see chapter 3, see final entry in Table 3-3). 

When the fuel flow rate smaller than 2.88 x 10° eyomm sec, a 
steady state diffusion flame is not obtained. This calculation 
was performed for Equation (3-9) without the d/dt term. The 
boundary conditions are given by Equations (3-11) and (3-12). 
Both 6 and en are functions of the fuel flow rate; for higher 
fuel flow rates , 6. is greater than the weap. When the fuel flow 
rate is smaller than 3.127 x 1094 cae sec, the reduced Damkohler 


number 3 becomes smaller than the extinction Damkohler pe 


indicating extinction. 


From the above analysis, it seems reasonable to say that for 


sustained pilot ignition, the Damkohler number for the corresponding 


boundary condition should be larger than the extinction Damkohler 


number, 


TEfSUEKSS 


which gives a conservative estimate for the minimum fuel flow 


Chapter 5 


THEORETICAL MODEL FOR PILOTED IGNITION OF WOOD. 


5.1 Background 


Ignition refers to the appearance of a flame in the volatile gas 
stream evolved from a solid exposed to external heating (usually 
radiative). Depending upon whether the ignition occurs with or without 
the aid of an external ignition source, it is accordingly classified as 


spontaneous (auto-) or piloted (forced). 


Several ignition criteria have been proposed, such as critical 
surface temperature at ignition [Simms (1963)], critical mass flux 
[Bamford et.al. (1946)], critical char depth [Sauer (1956)], critical 
mean solid temperature [Martin(1965)], etc. Of these, critical fuel 
mass flux at ignition seems to be physically the most correct since 
it can be related to flammability limits. However, surface temperature 
has proved to be the most useful ignition criterion since it can be 


conveniently related to fire spread [Atreya (1983)]. 


The ignition of cellulosic materials is a complex process 
involving both the gas and solid phases. When a cellulosic material 


(such as wood) is exposed to an intense radiant heat flux, the solid 
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chemically decomposes to inject fuel gases into the surrounding air 
and produce a flammable mixture which can be ignited by an ignition 


source. Therefore, this process involves flames of both pre-mixed and 


diffusion type. 


Kashiwagi (1974) suggests a theoretical model for auto-ignition 

based on the assumption that ignition occurs when the reaction rate in 
_5 3 

the boundary layer exceeds a value of about 10 g/cm sec. Gandhi (1986) 
suggests another model where it is postulated eee the auto-ignition 
point is the gradient reversal in the gas temperature profile at the 
solid-gas interface. Amos and Fernandez-Pello. (1988) suggest a 
theoretical model for auto-ignition by considering the absorption of 
radiation by the fuel vapor as a potential atte of ignition ss iiews 
model also considers the more practical case of the existence of 
convective flow over the combusting surface. Their analysis of this 
evolution process provides information about how a diffusion flame is 
generated from a combustion reaction which has initially premixed 
character. A similarly detailed analysis of piloted ignition has not 
appeared in the literature. Recently, Tzeng et.al.(1990) have conducted 
an investigation of piloted ignition. This predominantly gas-phase 
investigation did not consider the transient solid-phase decomposition. 


However, it showed the development of a diffusion flame in an initially 


premixed gas. This analysis is presented in chapter 3. 


The decomposition products of wood are very complex. Schwenker 
and Beck(1963) have detected 37 volatile compounds, Goos(1952) has 


identified 213 compounds as products of wood decomposition. The 
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composition of the pyrolysis nuddnere is also a function of the extent 
of decomposition, the initial density of wood, the temperature at which 
decomposition occurs, the fraction of char formed, the oxygen 
concentration in the air, the ambient pressure and the sample moisture 
content. Abu-Zaid (1988) and Nurbakhsh (1989) investigated the 
composition of pyrolysis products of wood. They found that the 
combustible portion is only about 25% (mass fraction) of the total 


pyrolyzed gas. 


5.2 Mathematical Model for Piloted Ignition of Cellulosic Solids 


Consider a thick slab of wood placed in an air stream under 
ambient conditions (as shown in Figure 5-1). The bottom face is 
considered impervious to both heat and mass transfer and the front 
face is exposed to a constant radiant heat flux. A part of this 
incident radiant energy is absorbed at the surface. As the solid 
surface temperature rises, it radiates a fraction of the energy back to 
the surroundings. Another fraction is convected to the adjacent air. 
The rest of the energy is conducted into the solid. Upon further 
heating, the solid undergoes thermal decomposition. The products of 
decomposition escape from within the solid through the surface into the 
gas. These pyrolyzates mix with the heated surrounding air in the 
boundary layer. Conditions at a distance h above the slab are 
maintained constant by a flowing oxidizer stream. A plane ignition 
source is placed near the porous plate. This ignition source is 


periodically "turned on" to test for piloted ignition. Mathematically 
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this process is described by the following equations. 


GOVERNING EQUATIONS: 

Gas Phase Equations. 

The chemical reaction is assumed to be a simple second-order, one- 
step irreversible Arrhenius reaction, Fuel + vO, + Products + q (heat) 
The governing equations are 

Energy equation: 

oT aT | a AnOma, Or 


Geert Vinx hoe SG ee acti (5-1) 
Bet Ou oex Cds CG 
P P 


Conservation of species 


Fuel: 
Oud a Y-¢ rs) a Ye 
Pp + pv came PD, J ee, (5-2) 
ot JFx ox 6.x 
Oxygen equation: 
ony: ony a a Yo 
Pp + pv =ascir (ed, ene (343) 
oat raha dx aix 
where: 
2 
w = Ap Y,Y EXP[-E/RT], 
sae at (5-4) 


The initial and boundary conditions are: 


Someta OOO << he 
Coron) 
See SY ome Yo i) and) eat Lane 
i ce) on foe) 
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at t = 0, the radiation source is turned on, resulting in 


production of the fuel. Thus, the boundary conditions after t = 0 


become: 
At xfs) 0: 
a Yp 
pd, = pv(Y, ¥ Ypo)> 
Os 
aay | (5-6) 
pD wl PMN Ec 
ra ibe : 
He 
S 


where pv, Yeu: 


the solid phase. 


ACRE xXee! Diet nO. 


(397) 
Ape CES OE yo ae | Ns 
12 re) oo Pr) 
Solid Phase Equations: 
Mass balance: 
OM a?) 
Wa Ww (5-8) 
ax Onc , 
Energy balance: 
rele gy 3 a Le 
[see G on aN ) (3-93 


and T. are determined by the heat and mass balance of 


ral 
Decomposition kinetics: 


dp 
WwW 
Pea Sesh, ~)p-1.)) PAPC Shean (5-10) 


The boundary conditions and initial conditions are: 


T(x, 0) = T (oe, t) ithe 


(Jari) 
pee) mip, MC, tee 
and the net heat flux into the solid is given by 
@ T(0,t) h AMS 
oc ek SIS a[ake- —(T, seers, meer o(T. raking} (5412) 
‘ee Bes € 


The gas phase equations are simplified by introducing the Howarth 


ak x 
BEans otmation,..£.=°——— p dx , t’ = t. They are nondimensionalized 
Poh, 0 


prAt Po Vv a 
pvececining + =——_z—;, M 6= , where T, is the 
Ee h pr T,-T 


fe) 
adiabatic flame temperature. Also, Be = E/RT,, Coe 1-(TL/T,), 
prp°ka°, Q-q/C (T,-T_), D= AC p_h_ EXP(=p°)/) 
? 8) £ © ? 2) fe) fe) ? 
PAGE i) 


LO EXP | 6 ara Finally, after assuming a Lewis 


and R = DY . 
1 - aid - 6) 


2 
number of unity and p D = constant, the energy equation and species 


conservation become 


6 
Li Y = wk’ Sous toll 6 (2-13) 
> ¢ 


Th. 


Meo Niko) © aires (5-14) 
s Fi 2 
OT rok 0€ 


where L(*) = 


The boundary conditions are: 


At fies ri Sele 
Oye B., 
s 
dY ,/é é= M(Y, - Yeo)? (5-15) 


dy (/a €=M Yo: 
andeat mec = mre 


Ye =20- (5-16) 


where M is determined by the mass balance at the solid gas interface. 
The corresponding non-dimensional solid-phase equations are presented 


in the next section. 


Gas-Phase Parameters for the Model 


In Chapter 3, the properties of the gas phase were assumed to be 
the same as that of a methane-air mixture. The activation energy for the 
gas phase reaction was assumed to be 29.1 Kcal/mole (as determined by 
Coffee, 1983). While this activation energy yielded good qualitative 
results, the results did not compare well Aerials eat with the 


ignition data for wood. Recently, Puri and Seshadri (1986) have studied 
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the extinction of a methane air diffusion flame. The activation energy 
calculated from the experimental data is around 45~47 Kcal/mole. Also, 
Westbrook and Dryer (1981) found that both higher (48 Kcal/mole) and 
lower values (30 Kcal/mole) of activation energies can predict the 
premixed flame velocity quite satisfactorily if the pre-exponential 
factor is adjusted. Now, since the flame velocity in a stoichiometric 
methane air mixture is 38 cm/sec, the pre-exponential factor is adjusted 
to pA = OMOEETION (cee yr The rest of the parameters (pi, are Ga Ty, T.) 
are the same as before (Chapter 3). (pr = 5.61x107- eal ais sec, 


Gs = 0.323cal/gm K, q = 11355 cal/gm fuel, T, = 2230 K, Tats 298K). 


The effect of higher and lower activation energies on piloted 


ignition has been investigated. The result is shown in Table 5-1. 
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Table 5-1. Effect of the higher (48 Kcal/mole) and lower (30 
Kceal/mole) activation energies. 


(Assuming the stoichimetric flame velocity is 38 cm/sec) 


Min. fuel flow 
rate for pilot 
ignition 


mS 
IS 


Flame temperature 
at min. fuel 
flow rate 


Activation 
energy 


Extensive studies on extinction of opposed-flow diffusion 
flames by Ishizuka and Tsuji (1981) have shown that near extinction, 
the limit flame temperature is remarkably constant for most hydrocarbons 
and is about 1550 K. Thus, it seems that the higher activation energy 
for methane is closer to the experimental data for diffusion flames. 
Since piloted ignition essentially represents the minimum requirement 
for the existence of a diffusion flame, the higher activation energy is 


used for the rest of this study. 


With these parameters, the above partial differential equations 
and the boundary conditions are solved by a implicit finite-difference 


method. 


ifs. 


5.3 Analytical Solution for Pyrolysis of Wood 


An analytical solution of the solid-phase temperature field, 
heat and mass transfer during piloted ignition of cellulosic solids was 
derived by Atreya (1983). This solution is used here and is outlined 


below. 


Consider a semi-infinite solid initially at constant ambient 
temperature LG At t > 0, the solid is exposed to a constant heat 
flux, F. Surface oxidation and internal heat transfer between the 
pyrolysis gases and the solid are ignored. The thermal properties are 
assumed to be temperature-independent and the net heat required to 
thermally decompose a unit mass of wood is taken to be zero. Since the 
heat transfer through an inert solid is by conduction only, the energy 


balance is described by: 


06. ) A 
ee ——— oie fOr xt sUemt > 0, 
wow wW 
ot Ox 
where a = Ty whee: T is the temperature of wood. The solid is assumed 


to be atthe temperature of its surroundings prior to the experiment. 
Therefore, the initial condition is @ = 0 for t < 0, x > 0. The boundary 
condition is obtained from the energy balance at the solid surface and 
is given by 

06 


-\ — = £ (6 Pea ES Ss ph se 
Ox x=0 
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where € = the emissivity (or absorptivity) of the surface and 
h = heat transfer coefficient/e, o = Stefan - Boltzmann constant, 


F 


the prescribed incident heat flux, sh = T. - To where T. is the 
surface temperature, f. = heat flux going into the semi-infinite solid 
at the surface, and Pras CA, are respectively the density, heat 
capacity and thermal conductivity of the solid. Because of the highly 
non-linear nature of the problem an exact analytical solution cannot 
likely be found. An approximate analytical solution for the above 
problem has been obtained by Atreya (1983) by the use of integral 


methods. Without entering into the details [see Atreya (1983)], the time 


as a function of surface temperature is obtained as 


2 Es 
6 r (250 +r-/B) (r+/By °° 6 (r¥256) 
s s s s 
a i aoc [a Lal, Nena ee ohaeebeea \ aedheaoalP ern 2) 
zon B (280 _+r+/B) (x-J/B) BE. 
2a pace _ 5 25 > "7 
where K = -— — r=- (h+ ase Des tp s=-— oT. , pp = ees). 
5 € 3 
* f. 
and fe =—=F+ré + s6 
s s s 


€ 

In the calculation of the surface temperature, a correction for 
the energy required to evaporate moisture was made. This quantity was 
estimated by using the constant weight-loss rate and latent heat 
of evaporation for water. The result was then substracted from the 
incident heat flux to obtained the corrected heat flux. This procedure 
is mathematically correct. The comparison between the experimental 
surface temperature and the calculated temperature as shown in Figure 


5-2, however, is better than expected. 
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The decomposition of wood is assumed to occur according to a 
first-order global chemical reaction with a known and fixed char density 


the mass balance of the one-dimensional decomposition process is 


Pr 

given by 
OM F op. 
Ox ot 


The decomposition kinetics are assumed to be described by 


dp 
Ww 
ae ey Le Re Pg) EXP (- E/T): 
at 
and the initial and boundary conditions are 
p(%,0) = poo. 
M (oe, t) = 0. 
To nondimensionalize the equations, a length scale is defined as 


L = X georeo/ &F which is obtained as the ratio of thermal conductivity and 


incident heat flux. Also, a diffusion time is defined as 7 


2 
1 [0X Po®) « The other non-dimensional variables are defined as follows: 


* * * * * 
T = T/T. Pp = PLS Poo x = x JL, t = t/t, E = EL /RT,: 
A ota e, S 

wd 


* ai fe) ¥ 


4 
6 a Piel Poo! d= ae yh) Ge M ty Pol oes oT. pa es 


* * 
r, = sed Aas 6 = alba 
Hence, the nondimensional mass balance and the decomposition kinetics 


equations are given by 


* * 
oM dp 
= - 
Ox 6t 

* 

and op 


* * ae 
ae = -A (p - § ) EXP(-E yak ys 


de 


Upon integration, 
Loe) 
* * * we * 
M. mien CL = 6¢) EXP[-E (@ + 1)]dx , 
0 
* 
where M. is the mass flow rate at the solid surface. During the initial 
* 
decomposition process, the density is almost constant, and p is nearly 
unity. Let 


Piene(H 4-45). BR Osa) 


* *2 
and @ = 1 + AAG + BBé 
S s s 


* * * * 
Forasmall.x;r.«.6 «= Gs - x eo; therefore, the upper limit of integration 
* 


* 
is replaced by x = dol /®.. By substituting this into the above 


ax 


equation, integrating, and keeping only the lowest-order terms of the 


* 
asymptotic expansion for large E , one obtains [Atreya and Wichman 


(1989) ], 
* * * as 
: A (1-6 ) Bb; < & EXP[-E (1-1/T. )] 
M. = ear T. EXP(-E /T. dy} {1 - i A We (5-18) 
S Ss 


Solid Phase Parameters 


The solid phase parameters have been determined by the parameter 
estimation technique from the experimental data for piloted ignition of 
Abu-Zaid (1988). The calculations are presented in Appendix D. The 


physical parameters used for wood are listed below: 


= 2 
h (convective heat transfer coefficient) = 2.0 W/m Kk, 
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A_ (pre-exponential factor for wood) = Sorte see 

A. (thermal conductivity of wood) = 0.167 W/m K Ca 
p._. (density of the wood) = 0.54 BS ot sta Snip 

E_ (activation energy for wood) = 31 Kcal/mole, 

C_ (heat capacity of wood) = 0.33 cal/g K, 

6. (char yield) = 0.25, 

€ (surface emissivity) = 0.75, 


T (ambient temperature) = 298 K. 


5.4 Numerical Prediction of Piloted Ignition Utilizing Analytical Solid 


Pyrolysis Solution. 


The numerical solution procedure for the gas phase is the same 
as before (Chapter 3). The surface temperature of the solid wood is 
calculated by using equation (5-17), and the fuel evolution rate is 
calculated by using equation (5-18). It is tactitly assumed that the 
contribution of any gas-phase exothermic reactions to the rise in the 


solid surface temperature is negligible during piloted ignition. 


The uniform grid spacing used in this work is 0.055 mm. A smaller 
grid spacing was also tested; Figure 5-3 shows the comparison the 0.055, 
mm grid spacing with the 0.027 mm grid spacing. The heat flux is 
assumed to be 1.88 Ween and the fuel is assumed to be 100% methane. 
Clearly, the temperature fields near the surface for small and regular 
grid spacings are very similar. The time for start of flashing for the 


small grid is 64./ sec, while for the regular grid it is 61.7 sec sane 
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time for sustained piloted ignition for the small grid is 71.1 sec, 
while for the regular grid is 70.7 sec. It is noteworthy that the 
ignition source has been turned on every 0.5 sec; thus the error is 
within one time period of turning on the ignition source. It is 
concluded that the grid space used in this study is appropriate for 
piloted ignition. Furthermore, the most important aspect of piloted 
ignition is the development of a sustained diffusion flame after 
the pre-mixed flame. The smaller grid can only predict a more accurate 
pre-mixed flame velocity. It will not significantly affect the time for 


sustained piloted ignition but will dramatically increase the CPU time. 


Figure 5-4 Shows the temperature history of the gas 0.15mm above 
the surface. The temperature peaks are due to the flashes and their 
subsequent quenching. Although it is not immediately obvious from the 
graphs (due to scaling), it is found that for small heat fluxes, 
the flashing periods are longer, while for higher heat fluxes, the 


flashing periods are shorter. 


5.4.1 Dilution Study: 


As mentioned earlier, numerous products are formed during the 
decomposition of cellulosic materials. For temperatures between 280°C 
to 500 mee the primary pyrolysis products are water, CO., CO, 
hydrogen, methane, ethane, acetic and formic acids, ethanol, aldehydes, 
ketones, tar and char. Abu-Zaid (1988) have analyzed the time- 


integrated product mass and composition of pyrolysis products for 
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FIGURE 5-4 Temperature history of the gas near the surface. 


(Assume fuel is 100% pure ) 
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Douglas fir. He found that the total hydrocarbons and carbon monoxide 
are only about a quarter of the total gas flux from the solid surface. 
Thus in this dilution study, the fuel concentration of the pyrolysis 
products is assume to be equal to 25%. Figure 5-5 is the temperature of 
the gas near the solid surface for the dilution study. The behavior is 
basically the same as that for the pure fuel case, only the time for 
flashes and ignition is a little bit larger. Table 5-2 shows a 
comparison of the ignition delay times for the experimental, pure fuel 
assumption and dilute fuel assumption. From Table 5-2, it can be seen 
that the ignition delay for the pure fuel assumption is about 20~30% 
lower than the experimental value, while the dilute fuel assumption is 


very close to the experimental value. 


Table 5-3 is a comparison of the fuel evolution rate for the 
pure fuel assumption and the dilute fuel assumption. The experimental 
fuel evolution rate is difficult to obtain from the experimental 
measurment, owing to the weight loss measurement errors. Roughly, the 
experimental measurment of the fuel evolution rate at piloted ignition 
iSpabOuUceO> ZecOnU se) ne cube As can be seen, the fuel evolution rate 
for pure fuel assumptiom is too low, while the result for the dilute 


fuel assumption is closer to the experimental data. 


Table 5-4 shows the comparison of the surface temperature at 
ignition for experimental, purefuel assumption, and the dilute fuel 
See er The surface temperature at ignition for pure fuel assumption 
is lower than the experimental value. Although the surface temperature 


for the dilute fuel assumption is also slightly lower than the 


TEMPERATURE 


TEMPERATURE 


us 


0.2 


0.1 


0.0 
ie) 


85 


WJ 
[oa 
E 
oe 
WJ 
[at 
= 
lJ 
b= 
10 20 30 ) 10 20 30 
TIME (sec 
DILUTE PEL Bee) TION DILUTE FUEL Bec) PTION 
0.3 
ree Ke) 
e 
jem 
lJ 
a 
= 
FE o4 
— F=2.60W/em2 0.0 — Fe2.62W/em2 
10 20 30 40 so te) 10 20 30 40 r%0) 
TIME (sec ne TIME (sec 
DILUTE FUEL eS MPTION DILUTE FUEL ASSUMPTION 
0.3 
i) 02 
= 
< 
ce 
uJ 
a 
= 
a4 
0.0 — Fei B8W/cm2 
0 20 40 60 80 100 


TIME een 
DILUTE FUEL ASSUMPTION 


FIGURE 5-5 Temperature history of the gas near the surface. 
(Assume fuel mass fraction = 0.25) 
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experimental value, it is closer to the experimental value. Furthermore, 
the surface temperature for the dilute fuel assumption is within the 


range of experimental error. 


From Table 5-2 to Table 5-4, it can be seen that the ignition 
delay time, the ignition fuel evolution rate and the ignition surface 
temperature for the dilute fuel assumption (25% methane and 75% inert) 
is close to the experimental data. Although it is within the 
experimental errors, the experimental results are consistently under- 
estimated. Further dilution (20% methane and 80% inert; results listed 
in the last column of Tables 5-2 to 5-4), brings the calculations closer 
to the experimental values. This amount of dilution is within the error 


of the experimental measurements of Abu-Zaid (1988). 
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Table 5-2. Comparison of the ignition delay. 


External 


heat flux | Experimental] Pure fuel Dilute fuel 


Dilute fuel 
assumption assumption 


assumption 


(25% CH,) (20% CH,) 


Table 5-3. Comparison of the fuel evolution rate at ignition. 


External 


heat flux Pure fuel Dilute fuel Dilute fuel 


2 assumption assumption assumption 


W/cm (25% methane) (20% methane) 
2 


(mg/cm sec) (aevemieees (ne emteec) (mg/cm sec) 


(The experimental fuel evolution rate at ignition is 


2 
about 0.2 to 0.25 mg/cm _ sec) 
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Table 5-4. Comparison of the surface temperature at ignition. 


External 
heat flux Experimental| Pure fuel Dilute fuel Dilute fuel 
assumption assumption assumption 


a (25% methane ) (20% methane ) 


( C) ( C) ( C) 


5.4.2 Flashes and Quenching 


From the experimental observations, there are usually some 
flashes which are subquently quenched prior to sustained piloted 
ignition. The sharp rises in surface temperature prior to sustained 
ignition in Figure 5-6 are due to flashes. These flashes are also 
observed in the numerical simulation. The sharp rises in temperature 


(prior to ignition) in Figures 5-4 are also due to flashes. 


The time elapsed between flashes and sustained ignition 
diminishes as the heat flux is increased. Also, the time between flashes 
and sustained ignition is larger for the dilute fuel than for the pure 
fuel, i.e., the time interval over which flashes occur is smaller for 


the pure-fuel case. 
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5.4.3 Comparison with Experimental Data for Other Woods. 


Figure 5-7 shows a comparison of predicted and measured surface 
temperatures at ignition. The experimental data are obtained from Atreya 
(1983). The solid line is for the pure fuel assumption, while the dashed 
lines are for the dilute fuel assumption. From the figure it can be seen 
that calculations for both of the assumed dilution levels are within the 
bounds of the experimental measurements. Clearly, the pure fuel 
assumption under-predicts the ignition temperature. Figure 5-8 shows a 
comparison of the predicted and experimentally measured ignition delay 
times. Once again, the pure fuel assumption under-predicts the ignition 


delay time and the dilute fuel assumption is more accurate. 


5.4.4 Comparison of the Minimum Fuel Evolution Rate with the Fuel 
Evolution Rate at Extinction of a Steady State Diffusion Flame. 


The minimum fuel evolution rate of a steady diffusion flame 
(i.e., near extinction) is calculated with the assumption of constant 
piloted ignition solid surface temperature. The numerical solution is 
obtained by solving the steady state version of equations (5-13) through 
(5-16). The analytical solution is obtained by large activation energy 


asymptotics (Linan 1974). 


Tables 5-5 and 5-6 show the comparison of the fuel evolution 
rate obtained from numerical simulation of piloted ignition (transient 


calculation) with the analytical and numerical calculations of the 
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minimum fuel evolution rates at extinction of a steady diffusion flame. 
The diffusion flame is adjacent to a surface whose temperature is the 
nearly constant surface temperature at piloted ignition for the 
transient calculation. From these tables, it can be seen that the fuel 
evolution rate at piloted ignition is quite close to the minimum fuel 
evolution rate of the steady state diffusion flame (i.e., near 
extinction). Although the values are slightly larger there is a trend 
toward increase with increase in heat flux, at least for the values in 


the first column. 
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Table 5-5. Comparison of steady state minimum fuel evolution rate 
with the numerical simulation of the piloted ignition. 
(pure methane assumption) 


External Fuel flow rate at Min. fuel flow Min. fuel flow 

heat flux piloted ignition, rate (steady rate (steady- 
same as in Table state analytical state numerical 
5-3 calculation) calculation) 


2 2 "; 
mg/cm sec mg/cm sec mg/cm sec 


(The surface temperature for the steady state solution is the same 
as for piloted ignition) 


Table 5-6. Comparison of steady state minimum fuel flow rate with 
the numerical simulation of the piloted ignition. 
( 25% methane and 75% inert gas) 


External Fuel flow rate at | Min. fuel flow Min. fuel flow 

heat flux | piloted ignition, rate (steady rate (steady- 
same as in Table state analytical state numerical 
5-3 calculation) calculation) 


2 2 2 
mg/cm sec mg/cm sec mg/cm sec 


(The surface temperature for the steady state solution is the same 
as for piloted ignition) 
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5.4.5 Effect of Air Velocity and Oxygen Concentration 


The boundary layer thickness is directly related to the ambient 
air velocity. Because the theoretical model in this study is one- 
dimensional, changes in the boundary layer thickness are the only method 
of simulating the changes due to ambient air velocity. Thus, an air 
velocity of 0.1 m/sec is assumed to correspond to a boundary layer 
thickness of 1.5 cm and a convective heat transfer coefficient of 

2 
2 W/m K. These values correspond to the experimental measurements. 
An air velocity of 1.0 m/sec will then correspond to a boundary layer 
thickness of 0.6 cm and a convective nee transfer coefficient of 


2 
4 W/m XK. 


Figure 5-9 shows the ignition delay time plotted against incident 
heat flux for two different air velocities and 0, concentrations. The 
ignition delay for 15% 0, is larger than that for 21% 0, and the 
ignition delay time for 1.0 m/sec air velocity is even larger than that 
for 15% O,. These results agree with the experimental measurements of 


Abu-Zaid (1988) qualitatively. 


5.4.6 The Effect of Moisture Content 


The effect of moisture content manifests itself in changing the 


physical properties of wood and diluting the products of pyrolysis. Thus, 


96 


O91 


*“suoTqerqUaoU0D ue3fxo 
pue SETIPOOTOA Ae AUSerSFFTp oF owyy AeTop uoyzaypusyz 6-S 


%7O PUD AVOOJOA IID JO 499449 
(98S) elu} UuOoIUb 


OI O8 OV 


JO%ZLT:S/wyo:Auq 
CO%G |S /w rO:alG 
CO%LT:S/WOQ’ | :ANQ 


FUNSI A 


+ 


Be 


<= 
N 


= 
ns 


(ZW9/M)XNIE JOEY JUSp!oU| 


97 
in this theoretical model, the effect of moisture is studied by 
changing the thermal conductivity, the heat capacity and the density of 
wood, The physical properties of wood as a function of moisture content 
are taken from Parker (1988); these are listed in Table 5-7. The fuel 
concentration as a function of moisture content is calculated by the 


following expression: 


p -ap 
dr char 
“= wiphy falbedee £4 ts Ske 
Beg) = Yen (dry) # (——S ih 
Py P char 
Here Py is the sample density when moisture is included; hence, Pel dey 
Figure 5-10 is the plot of ignition delay as a function of heat flux 


for various moisture contents; as expected, the higher the moisture 


content the larger is the ignition delay. 


Table 5-7. Thermal properties of wood as a function of moisture 


fuel 
concentration 


content after Parker (1988). 
C (heat 


Be se 
cBpacity) 
faisweiere | oe | ker ea fon 


5.4.7 Correlation of Results 


A(thermal 
conductivity) 


In Equation (5-17), if the surface temperature at ignition, (@_), 
is assumed to be constant, then the following relationship between 


ignition time and the incident heat flux (F) can be derived, 
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Hemeynverte ~ S'L, (5-19) 
where 
Zapper 
K ae i 2 . (5-20) 
SMe 
.9 § 
For constant K and L, the relationship between F and t ~ is 


linear, JKO is the slope of the line, and L is the intercept. 

Here L may be interpreted as the minimum heat flux for ignition, since 
it is the heat flux requlrcetene an infinite ignition delay time. 
Figures 5-11 and 5-12 show plots of square root of ignition time vs. 
heetyrlux. In both Figures 5-11 and 5-12/ jthe intercept is L»= anyon. 
while the experimental data of Abu-Zaid (1988) show that the intercept 
is about 1 Wow. The reason for this discrepancy is because the 
analytical solid pyrolysis solution (Equation 5-18) is obtained by 
assuming a constant density for wood during pyrolysis. For low heat 
fluxes, a lot of char is formed prior to piloted ignition, and the 
density near the surface is considerably lower. This means that the 
analytical solid phase pyrolysis equation over-predicts the mass flux 
for low heat fluxes, resulting in a smaller ignition delay time. This 
motivated us to find a complete numerical solution for wood 


decomposition. 


5.5. Numerical Solution for Decomposition of Wood 


The one-dimensional solid-phase mathematical model for the 
decomposition process studied here is given by equations (5-8) to 


(5-12). In these equations Cc. and Ne are functions of density. The 
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density of the solid in the pyrolyzing zone is assumed to change 
continuously from the initial density of wood; Piog? to the final 
density of char, Poe: At any instant, a partially pyrolyzed element 
of wood in this zone consists of char distributed throughout the 
unpyrolyzed active material. Since zero shrinkage is assumed, all 
densities are based on the original cola of the wood element; 
thus ee ots + Bre where p. and pe are the densities of the active 
wood material and char respectively. At t = 0, Poca sep and at 
t=, p = Poe: The densities as and Pp, can be described by the 


W 


following equations: 


e Pa” Pw 
(pay ake a! a 
Pad ~ Pwf 
=: Pad ; Pry 
PY = ici I) 
Pad Pwe 


while the thermal conductivity of the solid in the pyrolyzing zone is 


assumed to be given by 


where x and ar are the thermal conductivities of unpyrolyzed wood and 
char, respectively. Also, the specific heat of an element in the 


pyrolysis zone is assumed to be 


where, See and oy are the specific heats of unpyrolyzed wood and char, 
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respectivity. 


The differential equations (5-8) to (5-12) were solved by an 
implicit numerical scheme, and the numerical code was checked with the 
numerical program written by Atreya [1983]. The physical properties used 
for the numerical calculation are ae AGEGER ro = 0.07 W/m K 
SA a Ci. C ee 0.165 cal/g K. The remainder of the parameters are the 
same as before. Figure 5-13 shows the ignition delay time calculations 
using the numerical code for the pyrolysis of wood. A comparison with 
Figure 5-10 shows that for higher heat fluxes, the ignition delay time 
is not much different than the previous result (in which an analytical 
solution for the wood pyrolysis was used). However, at low heat fluxes, 
the ignition delay time for the numerical wood pyrolysis calculations 
is almost four times larger than the analytical wood pyrolysis 
calculations. This is because a substantial amount of char is formed 
over an extended exposure at low heat fluxes prior to ignition, as 
already discussed. Figure 5-14 is the plot of wood density at the time 
of piloted ignition. From the figure, it can be seen that the density 
near the surface for the low heat flux is much lower than for the the 
higher heat flux. 

_0 § 

Figure 5-15 is the plot of ( ignition time) -Vs incident heat 

flux. The slope of the line is different for different moisture 


2 
contents, and the intercept is about 0.9 W/cm . This agrees with the 


experimental results of Abu-Zaid (1988). 


In summary, the one-dimensional piloted ignition model presented 
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here captures all the essential features of piloted ignition of wood. 


5.6 Sensitivity Study 


In the above theoretical study of wood ignition, some of the 
physical properties were difficult to estimate accurately. Thus, the 
pertinent question is: how will the errors in the various parameters 
affect the ignition calculations. In what follows, the influence 
of solid-phase activation energy, frequency factor and thermal 


conductivity were investigated. 


5.6.1 Effect of Solid-Phase Activation Energy 


To understand the effect of activation energy on piloted 
ignition, three different values were used in the calculations. Figure 
5-16 shows the ignition delay time calculations for E28 Kceal/gmole, 
31 Kceal/gmole and 34 Kcal/gmole. It can be seen that a 10% change in 
the activation energy significantly changes the ignition delay time. 
For higher activation energy, the ignition delay time is much longer 
than that for the lower activation energy. This high sensitivity to 
activation energy also suggests that piloted ignition measurements may 
be a very accurate way of determining an overall activation energy for 
the one-step solid phase decomposition model. 


.o § 
Figure 5-1/7 is a plot of (ignition time) vs heat flux. It 
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can be seen that the intercept is 0.6 wien for EO = 28Kcal/gmole; For 
EY = 31 Kcal/gmole the intercept is 0.9 Genie and for ie 34 Keal/g 
mole, the intercept is 1.4 yi The slope of the three lines is 
therefore also changed. This is because the surface temperature at 
ignition (which affects both the heat loss and the slope) is changed 
with change in activation energy. The surface temperature for the three 
cases was 350, 410, and 470°C. Since the intercept of the line is 
interpreted as the minimum heat flux for ignition, the higher activation 
energy yields a higher minimum heat flux; a 10% increase in the 
activation energy results in a 50% increase in the minimum ignition heat 


flux. 


5.6.2 The Effect of Frequency Factor 


In the literature [Atreya (1983)], the frequency factor for the 
pyrolysis of wood ranges from 6.0x10. to rake In this study, three 
frequency factors A =6 0x10, 0.6x10. and 60x10. were investigated. 
Figure 5-18 shows the ignition delay for various frequency factors. 

As expected, higher frequency factors yield lower ignition delay times. 
0m 

Figure "5-19°is the plot) of (ignition time) vs heat fluxes 
trend and the slope of the separate lines is similar to the activation 
energy case. The intercept for different frequency factors is different. 
For AL = Oper Paces the minimum heat flux for ignition 1s about 
oe for AL = ee Oxi 0n/ sec the minimum heat flux is 0.9 Wen 


7 2 
while the minimum heat flux for AL = 60x10 /sec is 0.6 W/cm . In 
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other words, a factor of ten reduction in the frequency factor increases 
the minimum heat flux at ignition by 50%, whereas a 10% increase in the 
activation energy has a similar effect. This is caused by the linear 
dependence of the the reaction rate on Avs and the exponential 


dependence of the reaction rate on Es 


5.6.3 The Effect of Thermal Conductivity 


The literature values for the thermal conductivity of wood range 
from 0.134 W/mK to 0.263 W/mK [Atreya (1983)]. In this study, the 
thermal conductivity was assigned the values ro = 0.25 W/mK, 0.167 W/mK 
and 0.111 W/mK. Figure 5-20 shows the ignition delay time for these 
various cases. Although the ignition delay time for lower heat fluxes 
is very large, the three curves tends to practically the same asymptotic 
value. Figure 5-21 is the plot of (ignition cine ve heat flux, For 
A 70.25 W/mK, the intercept is about 0.8 Woe for A470. 167 W/mK, 
the intercept is about 0.9 dypare and for ae = 0.111 W/mK, the 
intercept is about 1.0 apy The thermal conductivity is related to 
the slope of the line. From equation (5-19), for higher thermal 
conductivity the slope is smaller and the intercept [or L in equation 
(5-19)] will be lower. This is confirmed by Figure 5-21. Clearly, the 
higher the thermal conductivity, the larger is the ignition delay time 
and the smaller is the minimum heat flux for ignition. Physically, this 
means that as the thermal conductivity increases, more heat have been 


conducted into the deeper area, the surface temperature and the 


pyrolysis product decreases. 


Lio 2 

From Figures 5-17 and 5-21, it can be seen that for higher 
activation energy the ignition delay is longer, and the intercept of the 
line is also larger, whereas, for higher thermal conductivity, the 
ignition delay is larger but the intercept of the line is lower. This is 
because the ignition temperature for higher thermal conductivity is 
lower (for 470.25 W/mK, T ignition is about 400°C, for r70.167 W/mK , 
T ignition is about 410°C, and for A,70 111 W/mK, T ignition is about 
420°C). However, the ignition temperature for higher activation energy 
istlargem.  Cror E34 Keal/gmole T ignition is about 470°C, for E31 
Keal/gmole T ignition is about irc ae ps E728 Keal/gmole, T ignition 
is. about 350°C). These changes in the ignition temperature in 
conjunction with Equation 5-19 explain the observed trends of ignition 
delay for different activation energies and different thermal 


conductivities. 
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Chapter 6 


Conclusions 


This chapter presents a summary of the important results for 
different components of this thesis. Section 6.1 summarizes the 
conclusions of the analytical-numerical method for the chemical reacting 
flow, section 6.2 summarizes the conclusions of the analytical solution 
for the extinction of the steady-state diffusion flame, section 6.3 
summarizes the conclusions of the one-dimentional model of piloted 
ignition, and section 6.4 summarizes the conclusions of the numerical 
simulation of piloted ignition of wood. Some recommedations for future 


work are presented in section 6.5. 


6.1. Analytical-Numerical Method for the Chemically Reacting Flow 


nll teThe required CPU time for the analytical-numerical method is 
less than for the implicit scheme. The explicit scheme for 
the chemical reaction scheme is fastest, but it is less 


stable when larger time steps are used. 


6.1.2. The analytical-numerical scheme is very efficient for the 
piloted ignition problem. Piloted ignition involves both 
pre-mixed and diffusion flames. The pre-mixed flame 


requires very small time steps, while the time step for 
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117 
diffusion flame can be 0(100) times larger than for the pre= 
mixed flame. The implicit and explicit methods overflow very 
easily during the transition from pre-mixed flame to the 
diffusion flame, while the analytical-numerical method does 


not overflow. 


6.2 Analytical Solution for the Extinction of the Steady-State Diffusion 


Flame. 


The extinction limit for a steady diffusion flame can be 
derived by large activation energy asymptotics (AEA), which provides 
expressions for the ignition and extinction Damkoher number. By using 
a higher mass flux, the flame temperature and the flame Damkohler number 
are obtained; then by reducing the fuel mass evolution rate, the 
Damkohler number is again calculated. If the Damkolher number is smaller 
wee the extinction Damkolher number, no solution for the diffusion 
flame exists, and the flame will be extinguished. The minimum fuel 
evolution rate for piloted ignition is very close to the analytical 


solution for the steady diffusion flame extinction limit (+103). 


6.3. One-dimensional Model of Pilot Ignition. 


6.3.1 The simple one-dimensional numerical model of piloted ignition 
describes the basic strunture of the piloted ignition process, 


which includes a pre-mixed flame, the quenching of a flash, 
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and the development of a diffusion flame after the flash. 


6.3.2 The minimum fuel evolution rate in itself is not sufficient to 


predict the onset of sustained piloted ignition, and the heat 


losses to the solid surface play an important role. 


6.3.3 The optimum location of the ignition source is found to be the 


eventual location of the steady diffusion flame. 


6.4 Theoretical Model of Pilot Ignition of Wood. 


oe 


6 


ei ye 


. Although the thermal conductivity of wood is a function of 


the exposed heat flux, temperature, moisture content and 
density, the overall thermal conductivity can be determined 
by the least square fitting of the experimental surface 
temperature. The curve fitting uses the analytical solution 
for the surface temperature of wood derived by Atreya (1983); 
the calculated thermal conductivity of dry Douglas fir is 


within the range of values listed in the literature. 


The integral analytical solution for the fuel mass production 
rate in the initial stages of Pohonpdsd ict can be used to 
determine the pre-exponential factor of the pyrolysis of 
wood by using the least square fitting of experimental 
measured mass flux around piloted ignition. The predicted 


pre-exponential factor for dry Douglas fir is within the 
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range of values’ listed in the literature* ‘Since’ the 
analytical“solution is only’valid for the initial”séageuge 
decomposition, this model cannot be used for lower heat flux 


where ignition delay is very large and a lot of char has been 


formed before piloted ignition. 


By using the thermal conductivity obtained by the least 
square fitting, the analytical solution can predict the 


experimental surface temperature of wood very accurately. 


This numerical model for the piloted ignition of wood can 
predict the flashes, the quenching and sustained piloted 
ignition. For lower heat flux, the period of flashes and 
quenching is longer, while for higher heat fluxes, the 


flashes and the quenching are not so obvious. 


During the flash-quenching period, the fuel evolution rate is 
much lower than the minimum fuel evolution rate of the steady 
state diffusion flame, whereas at sustained piloted ignition 
the fuel evolution rate is very close to the minimum fuel 
evolution rate. This leads us to believe that conditions for 
piloted ignition for wood are quite similar to those for 


extinction of a diffusion flame. 


The assumption that the fuel evolute from the solid wood 
is the same as pure methane underestimates the ignition 


delay, the ignition surface temperature, and the ignition 


ae 


eo. 
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mass evolution rate. The predicted ignition delay is about 
30% lower than the experimental value, the predicted ignition 
surface temperature is about 15% lower than the experimental 
ignition surface temperature, while the predicted mass 
evolution rate only about 1/5 to 1//7 of the experimentally 


observed value. 


The assumption that the fuel evoluted from the solid wood is 
25% methane and 75% inert gas (by mass) enables the 
experimental ignition delay, the ignition surface temperature 
and the frat tion mass evolution rate to be accurately 
predicted. From the experimental observations, the gas flow 
out of the solid in the early stages of pyrolysis is 
primarily water vapor, CO, CO, CH), and higher 

hydrocarbons; it is evident that the concentration of total 
hydrocarbons and CO is very low at piloted ignition. Both the 
surface temperature and fuel evolution rate are very 
important for the piloted ignition, but the dilusion of 


the fuel is also a very important factor. 


The study of the effect of 0, concentration and ambient air 
velocity shows that higher air velocities and lower 0, 
concentrations increases the ignition delay and increase 
the minimum heat flux for piloted ignition. This result is 
agreement with the experimental observation of Abu-zaid 


(1988). 


6.4.10 


dy 
The study of the effect of moisture content shows that a 
higher moisture content increases the ignition delay without 


altering the minimum heat flux for ignition. 


The model which use the analytical solution for the wood 
pyrolysis equations [derived by Atreya and Wichman (1989) ] 
predicts an ignition delay at high heat flux (greater than 
2.6 wens very close to that calculated from a model which 
solves the pyrolysis equations numerically. At lower heat 
fluxes (less than 2.0 wrens) the ignition delay predicted 
by the analytical solid solution model is only 1/3 of the 
value predicted by the numerical solid-phase model. The 
DLOC OL <6 Seareron eee vs heat flux shows that the 
numerical solid phase model obtains closer value of 
minimum heat flux for piloted ignition to the experimental 
observation of Abu-zaid (1988). The analytical solid phase 


model underestimates the minimum heat flux for piloted 


ignition. 


6.4.11 A sensitivity analysis shows that: 


A. The activition energy for thermal decompostition of wood 


is very important for the prediction of piloted ignition. 
An 0(10%) change in the activation energy makes an 0(100%) 
change in the ignition delay. The plot of (t ienitiontae 4 
vs heat flux shows that the curves for different activation 


energies have similar slopes, but that the minimum heat 


fluxes for piloted ignition are different. An 0(10%) 
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increase in the activation energy makes an 0(50%) increase 
in the minimum heat flux. 

B. The frequency factor is very important for the prediction 
of piloted ignition, but is less sensitive than the 
activation energy. An increase in the frequency factor 

SOS 
reduces the ignition delay. The plot of (t ignition) 
vs heat flux shows that the curves for different frequency 
factors have the same slope, but the minimum heat flux 
Lorapilot ignitionsis difference, (This trend is similar 
to the effect of activation energy). The effect of a 
factor of .10 increase in the frequency factor is similar 
to reducing the activation energy by 103%. 

C. The thermal conductivity is very important for the 
prediction of piloted ignition; the higher thermal 
conductivity, the larger the ignition delay. The effect of 
thermal conductivity on the ignition delay is different 

er 

from the activation energy. The plot of (t ignition) — 
vs heat flux shows that the lines with higher thermal 
conductivity have a lower slope, and that a higher thermal 


conductivity produces lower minimum heat fluxes for 


piloted ignition. 


6.4.12 A new definition of piloted ignition is proposed in this 
work. It is found that the development of a sustained 
diffusion flame after the pre-mixed flash successfully 
predicts the experimentally measured ignition delay, ignition 


surface temperature, ignition mass flux and minimum fuel flow 


23 
rate for piloted ignition. This model is helpful in 
clearifying important parameters for piloted ignition that 
are not easy observed from the experimental data. (such as 
quenching effect, location of ignition source, thermal 


conductivity of wood and activation energy of wood). 


6.5. Recommandation for the Future Work 


In the present model of piloted ignition, the chemical reaction 
in the gas phase is assumed to be a one-step irreversible chemical 
reaction. In some pratical experimental observations, the one-step 
chemical reaction is too simple and multistep chemical reactions may 
produce more accurate results. Besides, the composition pyrolysis 
products are very complex, a great deal of CO and water vapor are 
contained in the pyrolyzate; thus, a multistep chemical scheme would 
give a better representation for the effect of water vapor and CO 


during piloted ignition. 


For the time being, the piloted ignition model in the gas phase 
is one-dimensional. The actural gas phase in the boundary layer is two 
dimensional. Since piloted ignition involves pre-mixed flames, and since 
the pre-mixed flame thickness is very small (only about 0.1mm), the grid 
spacing should be smaller than the flame thickness. A two-dimensional 
model for the gas phase would be more practical, but such a model might 


require a more powerful tool, such as a super computer. 


Appendix A 


A COMBINED ANALYTICAL-NUMERICAL SOLUTION PROCEDURE FOR 


CHEMICALLY-REACTING FLOWS 


A.1. INTRODUCTION 


The possibility of exploiting the often vastly different time 
scales for chemical reaction [0(1Obeseddl and for Rede conte 
convection and diffusion of this chemical release into the flow field 
[O(1 sec)] appears to have been first discussed physically by Baum and 
Corley(1981). (A more detailed history of the "operator-splitting" 
method is given in the recent study of Wichman (1990)). In their short 
work Baum and Corley (1981) analyzed the laminar mixing of fuel and 
oxidizer streams in a. plane channel flow for arbitrary Lewis numbers, 
and for two reaction models; (1) delta-function heat release, and 
(2) one-step Arrhenius kinetics. For (2), the solution is obtained by 
dividing each streamwise marching step into two parts. In the first 
part, a locally homogeneous rate chemistry problem is solved by 
quadratures (for (1) the analytical solution of this part is trivial). 
This solution is then applied as a (local) initial condition for the 
second part, in which this heat release is transported by convection 


and diffusion to the neighboring fluid elements. Thus, the solution is 
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obtained by an "operator-splitting" procedure. 


The purpose of the present study is to develop a rational 
motivation for the use of the operator-splitting procedure in combustion 
theory by analyzing a simple model problem in some detail. A numerical 
scheme that is consistent with this model is then applied to a simple 
premixed laminar flame model that has been previously employed as a test 
for various numerical integration schemes [Peters (1982)], such as the 


method of lines, adaptive grids explicit and implicit schemes. 


A.2. The Model Problem 


Consider a purely thermal model, ton which the gas-phase energy 
balance is represented by transient, diffusion and heat-release terms. 
When properly nondimensionalized, the coefficients of the transient and 
diffusive terms are unity, while the heat-release term contains the 
factor ID = At ,/(4t.), called the Damkohler number. Here at, is the 
characteristic diffusion (or convection) time, while ae is the 
characteristic reaction time. Thus, we consider the following initial 


boundary-value problem, 


2 
oT oie 
Lim etd DOK AM -o <x <=, t>0 
ot Ox 


(A-1) 
Taye amet (op oy eas0s 


ll 


T(x, t)) TO (x). 


Note that the heat-release term (analogous to the chemical reaction term 
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in the combustion equations) is written as a simple function of x and t. 
The solution for T(x,At) can therefore be obtained from the solution at 
the previous time step ch from the equation (we put 45 = 0 with no loss 


of generality) 


At ja fore) 
T(x,At) = - G(x,At|é,r) [DQ(&,7r) ]dédr - TO (6) G(x,4,t|é,0)dé, (A-2) 
0) -<0 -2 
where 


1 cel Bee 
EXPY Se eehee ty (oo 


2/n(at -. 7) 4(At - 7) 


G(x,at|é,r7) Pike 


is the Green’s function for equation (A-1). The first term in 
equation(A-2) describes the influence of the inhomogeneous reaction 
term, while the second term describes the redistribution of energy from 
the previous time step by the homogeneous diffusion operator 


2 2 
d(e)/dt - d (e)/dx = 0. Note that the second term does not contain Q. 


The difficulty with obtaining numerical solutions of equations 
such as equation (A-1) arises from the largeness of D. This usually 
necessitates taking very small time steps, of order At. Rational means 
for avoiding this requirement (and therefore speeding up the 
computation substantially) are sought by analyzing equation (A-2) as 


D+ o@, 


In order that all three terms in equation (A-2) balance as D + o, 
the time interval At in the first integral must become very small; in 
eel 
fact, At must be of order D = At.. Thus, by defining s = Dr as a new 


0(1) variable of integration, one finds, as D + o~, 
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foe) foe) 


T(x,At) = - Dat G(x, At|é,0)Q(E,0) dé - T (€)G(x, at] é,0)dé, (A-4) 


-co -@ 


where At is now O(at.). For very small At the Green's functions in 


(A-4) reduce to the Dirac delta function -6(€-x), whereby 


T(x,At < At.) = at DQ(x,0) + T (x), (A-5) 


A 
> 
ct 
R 

S 


which is valid for At < 


For the remainder of the time step a finite solution for T is 
possible only if Q(x,t = At.) = Q. This means that equation (A-1) 
describes a transport-limited process (convection,or diffusion, or 
both), in which an increment of combustible gas is transported to 
the reaction zone, where it reacts "instantaneously", thereby producing 
heat and depleting the increment of combustible gas. The heat generated 
by reaction is transported away from the reaction zone, and further 
reaction cannot occur [i.e., Q is now zero] until another increment 
of combutible gas has been transported into the reaction zone. Thus, 
the overall time interval for this process is approximately one 
bLee Sot aeat 


characteristic transport time At in equation(A-2). ] 


d d 


Since chemical reaction and heat release occur at the beginning of 
this time interval, equation (A-5) may be used as an initial condition 
for the calculation of the change in T over the remainder of the time 


step. Thus, equation (A-2) gives 


T(xeAL) a= ae T(x, At )G(x,at|€,0)dé, (A-6) 
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where T(x, At) is given by equation (A-5), and G may be evaluated at 


Lewes inetead of at-at—<-nt 
c 


In sum, the solution of equation (A-1) over the time interval 
At can be obtained by first solving the spatially-homogeneous equation 
dT/dt = DQ, giving T(x, At.) = TS (x) + at. DQ(x,0)[which is equation 
(A-5)], and then using this as a modified initial condition for the 
solution over the remainder of the time step. Then one solves dT/dt = 


2 2 
6 T/dx , obtaining equation (A-6). 


The preceding physical and mathematical discussion implies that 
this integration procedure can be implemented wherever the local 
reaction rate quickly consumes the local excess combustible. If. this 
were not. the case, then the proposed integration scheme could be applied 
only to diffusion flames with infinitely fast chemistry in 
infinitesimally thin reaction zones. In general, computational methods 
are unnecessary for such problems. Consequently, this integration 
procedure is expected to apply for diffusion flames with "slow" 
chemistry and premixed flames; for diffusion flames the excess component 
is fuel or oxidizer, while for premixed flames it is the combustible 
mixture itself. Such a general level of applicability is important 
because most currently challenging problems in combustion theory involve 
flames of mixed character. Therefore, the first test of this solution 


procedure is a premixed-flame-speed calculation. 


A.3. The Premixed Laminar Flame 


The premixed flame equations solved here are 
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2 
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where 8 is the Zeldovich number, and D is the "Damkohler number". Over 
the reaction time interval one solves dT/dt = DR, dY/dt = -DR. The 
Schvab-Zeldovich procedure gives Y = C(x) - T, thereby requiring the 


solution of only one nonlinear equation, viz., 


aT A(1-T) 
— = D(C(x) - T)EXP}|- ——————— . (A-8) 
at 2 aChis. i 


This equation was integrated with a fourth-order Runge-Kutta scheme over 


: s sah te n++ n n+1 n ; 
the entire time step, At, giving T, = Teh ; Y; = Yon , which were 


then used as augmented (local) initial conditions for the final 
A n+} n+? : : : 
calculation of T, ; Y; from the homogeneous diffusion equations 


[put R = 0 in Eqs. (A-7)]. This calculation employed a simple (implicit) 


centered-difference scheme, viz., 


l 1 i 1 
i) T ee alT. gy _ oT n+ iT Og _T. n 
af i+ i- ih 
& 1 1 1 (A-9) 
ii) y mt cine pts _ oT n+ Days tech _y n 
i ma i+ ib i- ih 


where a = At/Ax 


The parameter values used were Le = C054 WSO!) DIG sel CLOSE es 
Ax = (0.25, 0.1, 0/0625), and At = (0700055,0.00), 0.002..0.005 =u aia 


2 
0.02, 0.05). Note that asymptotic analysis gives D = B /2Le; note also 
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2 
tnaGepevele -- O(1), as required for the anaylsis of Sec: I1. 


A.4. Benefit of the Analytical-Numerical Method 


The numerical integration scheme of Sec. III was compared to the 
flame-speed calculations presented in the GAMM workshop [Peters (1982) ], 
and to calculations of the same problem using an explicit numerical 

1 
scheme and an implicit numerical scheme (iterate with the reaction 
term). An example of these latter comparisons is shown in Table 1 for 
the case Le = 0.5, 6 = 10. The numerical scheme of Sec. III converges 
to a solution in all cases, whereas the other schemes develop 
convergence and overflow problems. The required CPU time is less than 
for the implicit scheme, though the explicit method is fastest when it 
works. In multiple space dimensions it is believed that his method will 
be much quicker and easier to implement than other, more sophisticated 


algorithms. In addition, flames of mixed (diffusion/premixed) character 


are not expected to produce excessive computational difficulties. 


Remarks: !. The scheme of Sec. III is essentially "adiabatic", while the 
explicit scheme is “constant temperature". Therefore, they 
represent, in a sense, thermodynamic extremes. Explicit 
schemes have also been used by Reitz(1981). 


fey 


Table A-1: Flame= Velocity for Le=0.5, #=10, V=flame velocity, 
N=number of time step calculated, (CPU time, by Digital 
MiOronVAna Ll w unit sec), 
t =, 0806t 0.002 0.005 0.01 O02 02805 Ov) 


method 1, case 1 
V=0.962 V=0.960 V=0.945 V=0.944 V=0.927 V=0.850 V=0.624 


N=900 N=450 N=150 N=100 N=120 N=40 N=40 
GPUCTAWe ROS 88.54 52 P17. pT aN Be 46.28 30a t9 70.54 


Ax=0.25 method 2, case 1 V=0.965 V=0.965 V=0.959 not 
N=900 N=450 #£xN=150 convergent 
GPU 204266 |17259). V62 216 


method 3, case 1 V=0.964 V=0.963 V=0.954 V=0.938 V=0.952 over 
N=900 N=450 N=150 N=240 N=50 
CPU Bo67, 21.69 LOSS 2 i363 6222 


method 1, case 2 
V=0.962 V=0.960 V=0.951 V=0.936 V=0.907 V=0.852 V=0.777 
N=900 N=450 N=150 N=60 N=30 N=40 N=40 
GPU 1425634 2140730 74074 324825 30256 87icBs L69.294 


Ax=0.10 method 2, case 2 V=0.965 V=0.966 V=0.968 not 
N=900 N=450 N=150 convergent 
CPU 58425 429.57 284.78 


method 3, case 2 V=0.964 V=0.964 V=0.961 V=0.955 V=0.945 over 
N=900 N=450 N=150 N=60 N=30 flow 
CPU 9224 48.32 eS a7, 9.95 6.69 


method 1, case 3 
V=0.962 V=0.959 V=0.952 V=0.936 V=0.906 V=0.853 V=0.777 
N=900 N=450 N=150 N=60 N=30 N=40 N=40 
CPU 678.98 BAZ Use, 2:)faSt 95 3 46°. 19 L138 4orerzou Le 


Ax=0.063} method 2, case 3 V=0.965 V=0.965 V=0.965 not 
N=900 N=450 #£xN=150 convergent 
GPU Tee 493 eae 2 oGaS 


method 3, case 3 V=0.963 V=0.963 V=0.960 V=0.954 V=0.944 over 
N=900 N=450 N=150 N=60 N=30 flow 
CPU LEO 72974 Zien pee a Onis 


V=0.944 (analytical solution) 

Method 1, analytical numerical method. 

Method 2, Implicit method. 

Method 3, Explicit the chemical reaction term and 
implicit for the convection diffusion term. 


Remarks: 
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APPENDIX B 


PROGRAM LISTING FOR THE PILOT IGNITION MODEL 


FH HK HAHAHAHA HAHN A HHA HA HHA HHA AHN A AERA HERE TE EK EH PER EATER TEE TER TEE TE TE Te aa GE ME 
COMPUTER MODEL FOR THE PILOT IGNITION OF WOOD 
by: Lin-Shyang Tzeng 
FIRDI in Taiwan 
Peroee box c4onnsinchu sUGs  aiwan Rh. OG: 
SEPTEMBER 1990 
HAHAHAHAHA AAA AHHH A HEH AHK HAHAH HHA HHA HHH HAKKAR ERATE RTE TE TE ME HE Mee 2 ak 
The following five lines is the example of input data (see next page 
under C “input data") 
peer > 0. 00002 7670,00027670,00666 0.0037 °0-00666 0.02 0 0 0 .218 1.0 1.657 
Com SILO ore obi mune mote OGL 98. 06205197 6.OR/ 75 0.232 
#eQ00 560 110 21 221 241 3 40,10 0.241 20 1 0 
set O97.0.69.0.0710 70183100 
PNM1F38 . DAT 
The case shown above is in the file named PNM1F38.DAT 
combined analytical & numerical method 
solve chem reac. analytically when T(I). gt. 0.95 
including the surface temp. 
solve the solid equation numerically 


not dumping the data after turn on the ignition source 
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ios 

Le number not = 1. 

PROGRAM ONEDPISN 

COMMON T(401) , YF(401) ,YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, NDUM 
1,NIGS,DTOD, IM, IN, 10, NODUF, TIME, DT, IFLAG, IDUMP,HL, NMAX, SMA 
1,1,DT1,DI2, DIB 9CE.CO ,ARSBT.O , DAGFS.REAGHCEAS CoMeSTY 
1,AFU(6),DF(6) ,BF(6) , YFS 
COMMON/WDATA/WTD(451) , WIDO(451) ,WAMD(451) , WAMDO(451) , WDN(451) 
1 ,WDNA(451),WDNC(451) ,WHCP(451) ,WCND(451) ,WA(451) ,WD(451) ,SIGM 
1 ,WD1(451) ,WDNO(451) ,WB(451) ,WR(451) , DNWF, DNF ,HCPA, HCPC, CNDC 
1 ,WDX,WTIME, EACT 

DIMENSION TOD(401) , YFOD(401) , YOD(401) ,TTD(401) , YFTD(401) 
1,YOTD(401) 

DOUBLE PRECISION WDN,WDNO,WDX 

CHARACTER*16 SFTPDAT 

OPEN (UNIT=10 , FILE=' IGDUFLN. DAT’ , TYPE='NEW’ ) 

OPEN (UNIT=11, FILE=’ DAISN.DAT’ , TYPE=’OLD’ ) 

OPEN (UNIT=12, FILE='OUPTN.DAT’ , TYPE=’ NEW’ ) 

OPEN (UNIT=13 , FILE=’HTFXN.DAT’ , TYPE=' NEW’ ) 

IGDUFLN.DAT data for check if there is ignition 

DAISN.DAT input data file 

OUPTN.DAT file for temperature and species profile 

HTFXN.DAT file for wave of heat flux 

input data 

READ (11, *) DT), DEZ#DT3) DXIS DK? DXB SSMi SMP Si ST) wr Seance 
READ(11,*)AF,BT,Q,DA,FS,HSF,CVF,TAMB,DNS ,RKK,APS,EPSL, YOS 
READ(11,*) NMAX,NIG, IG, IM, IN, 10, IGN, NDUP ,NHS , NW, IDUMP ,NCH.M 


READ(11,*)DNF,HCPA, HCPC, CNDC,WDX, EACT 
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READ (11,2400) SFTPDAT 
2400 FORMAT (A16) 
OPEN(UNIT=14 , FILE=SFTPDAT, TYPE=' NEW’ ) 
WRITE(14,*)DT,DT2,DT3 ,DX1, DX2 ,DK35SMD,SM2 ,ST1, ST2;YFS,, RLE 
WRITE(14,*)AF,BT,Q,DA,FS,HSF,CVF,TAMB,DNS ,RKK,APS,EPSL, YOS 
WRITE(14,*)NMAX,NIG,IG,IM,IN,IO, IGN,NDUP,NHS ,NW, IDUMP,NCH,M 
WRITE(14,*)DNF,HCPA, HCPC, CNDC,WDX, EACT 
WRITE(14, 2500) SFTPDAT 
2500 FORMAT(1X,A16) 
DT = time step of no chem. react. 
DT2 = time step of premixed flame 
DT3 = time step of diffusion flame 
DX1 = grid space near the surface 
DX2 = grid space in the middle area 
DX3 = grid space near the upper surface 
SM1 = const fuel flow rate 
SM2 = coeff. of linear fuel flow rate 
STl = constant surface temp. 
ST2 = coeff. of linear surface temp. 
SM1,SM2,ST1,ST2 in this program are only dummy variable 


YFS = fuel concentration in the solid phase 


RLE Lewis number 

AF = non-dimensional flame temperature 

BT = non-dimensional activation energy of fuel 
Q = non-dimensional heat of reaction of fuel 


DA Damkohler number 


FS stoichiometry ratio 


je 
HSF = radiation heat flux imposed on the wood 
CVF = convection heat transfer coeff. (W/m2 k) 


TAMB = ambient temperature k 


DNS = density of wood g/cm3 
RKK = thermal conductivity of wood W/m k 
APS = pre-exponential factor of pyrolysis of wood 


EPSL = surface emissivity 

YOS = oxygen concentration of ambient air 

NMAX = maximum number of step 

NIG = number of time step to turn on the ignition source 

IG = grid location of the ignition source 

IM = number of the grid near the surace 

IN = number of the grid in the middle area 

- IO = total number of the grid 

IGN = number of eride for ignition source 

NDUP = number of step between calling subroutine DUMP 

NHS = number of step between turn on the ignition source 

NW = flag to control if write down data of temp. fuel and time 
IDUMP = number of data grids to be dumped 

NCH = interval of step to write down the surface temp. with respect 
to time. 

M = number of step to calculate the chemical reaction term 


DNF = final char density 


HCPA specific heat of active wood 
HCPC = specific heat of char 


CNDC 


thermal conductivity of char 


WDX = grid space of wood 
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EACT = activation energy of wood (kcal/g mole) 
DX1* (IM-1)+DX2*( IN-IM)+DX3*(IO-IN) = 1.00 for boundary layer thickness 
= 1.5 cm 

TRES=0.0011 

TMAX=(DT+DT3) *0O. 5*NMAX+1.5 

initial condition 

DO 100 I=1,10 

T(1I)=0. 

YF(I)=0. 

YO(I)=YOS 

100 CONTINUE 

coefficient of matrix and boundary cond. 
B(1)=0. 
A(1)=(-1.)/DX1 
AFU(1)=-1./(DX1*RLE) 

NIGA=0 

NIGB=0 

NIGS=0 

NDUM=0 

NODUF=0 

NTD=0 

NTD1=0 

MA=0 

NIGT=NIG 

DTTD=DT 

NDUPTD=NDUP 


NHSTD=NHS 


an 
LIGN=IG+(IGN/2) 
TIGPS=0 
DTOD=DT 
TIME=DT 
ERRO1=0.0 


ERRO=1.E-8 


INOWR=0 
specify the x-coordinate XX(I) (Variable grid) 
XX(1)=0. 
DO 102 I=2,1M 
XX(1)=XX(1-1)+DX1 
102 CONTINUE 
DO 104 I=IM+1,IN 
XX(1)=XX(1-1)+DX2 
104 CONTINUE 
DO 106 I=IN+1,10 
XX(1)=XX(1I-1)+DX3 
106 CONTINUE 
send data to temporary file 
110 CONTINUE 
DODLO Se l= 110 
TOD(1)=T(1) 
YFOD(1)=YF(1) 
YOD(1I)=YO(I) 


108 CONTINUE 
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start the new time step 
111 CONTINUE 
IF(TIME.GT.TMAX) GO TO 400 
IF((MA-M).GT.2) GO TO 150 
NCH1=(N/NCH) *NCH*NW 
NCH2=(N/NCH) *NCH 
IF(N.EQ.NCH1) WRITE(13 ,1100)TIME,T(2) ,T(IG) 
eC), Yee) LE CLGyY VECrNy sce) 0 Le), YOCIN) 
CALL SFIN(SMA,STP,DT,TIME,HSF,CVF,TAMB,DNS ,RKK,APS , EPSL,N) 
IF(N.EQ.NCH2) THEN 
locate the max temp 
IDUF=2 
TDUF=0. 
DO 40 I=2,10-1 
PrCTCL Greer vury GO TO*S0 
GO TO 40 
30 CONTINUE 
IDUF=I1 
TDUF=T (I) 
40 CONTINUE 
TISEC=TIME/0.55187 
WRITE(14 , 2200) TISEC,SMA,T(4) ,T(1),TDUF, IDUF,N 
END IF 
2200 FORMAT(CIX ,“tsecyoMA Size" TPM Na FOS GEOS, LX Tas TX 15) 
IF(SMA.GT.2.) GO TO 400 
T(1)=STP 


YF(1)=SMA*YFS 


YO(1)=0. 

D(1)=SMA+(1./DX1) 

A(2)=DT*((SMA/DX1) - (1. /(DX1**2) )) 

D(2)=1. - (SMA*DT/DX1)+(2.*DT/(DX1**2) ) 
B(2)=(-1.)*DT/(DX1**2) 
A(3)=SMA*DT/DX2-2.*DT/((DX1+DX2) *DX2) 
D(3)=1.+2.*DT/((DX1+DX2) *DX2)+2.*DT/( (DX1+DX2) *DX1) 
1-SMA*DT/DX2 

B(3)=(-2.)*DT/( (DX1+DX2) *DX1) 

A(4)=(SMA*DT/DX2) -DT/(DX2**2) 

D(4)=1. - (SMA*DT/DX2)+(2.*DT/(DX2**2) ) 
B(4)=(-1.)*DT/(DX2**2) 
A(5)=SMA*DT/DX3-2.*DT/((DX2+DX3) *DX3) 
D(5)=1.+2.*DT/((DX2+DX3) *DX3)+2 .*DT/( (DX2+DX3) *DX2) 
1-SMA*DT/DX3 

B(5)=(-2.)*DT/((DX2+DX3) *DX2) 
A(6)=(SMA*DT/DX3) -DT/ (DX3**2) 

D(6)=1. - (SMA*DT/DX3)+(2.*DT/(DX3**2) ) 
B(6)=(-1.)*DT/(DX3**2) 

DF(1)=SMA+(1./(DX1*RLE) ) 
AFU(2)=DT*((SMA/DX1) - (1. /( (DX1**2) *RLE) ) ) 

DF(2)=1. - (SMA*DT/DX1)+(2.*DT/((DX1**2) *RLE) ) 
BF(2)=-1.*DT/( (DX1**2) *RLE) 
AFU(3)=SMA*DT/DX2-2.*DT/(( (DX1+DX2) *DX2) *RLE) 
DF(3)=1.+2.*DT/(( (DX1+DX2)*DX2)*RLE) +2. *DT/(( (DX1+DX2) *DX1)* 
1 RLE) -SMA*DT/DX2 


BF(3)=-2.*DT/(( (DX1+DX2) *DX1) *RLE) 


eZ 


140 
AFU(4)=(SMA*DT/DX2) -DT/( (DX2**2) *RLE) 
DF(4)=1.- (SMA*DT/DX2)+(2.*DT/( (DX2**2) *RLE) ) 
BF(4)=-1.*DT/( (DX2**2) *RLE) 
AFU(5)=SMA*DT/DX3-2.*DT/( ( (DX2+DX3 ) *DX3 ) *RLE) 
DF(5)=1.+2.*DT/( ( (DX2+DX3) *DX3) *RLE) +2 .*DT/( ( (DX2+DX3 ) 

1 *DX2)*RLE) -SMA*DT/DX3 
BF(5)=-2.*DT/(( (DX2+DX3) *DX2 ) *RLE) 
AFU(6)=SMA*DT /DX3-DT/((DX3**2) *RLE) 
DF(6)=1. - (SMA*DT/DX3)+(2.*DT/( (DX3**2) *RLE) ) 
BF(6)=-1.*DT/( (DX3**2) *RLE) 
calculate the homogeneous chemical reaction 
DT1=DT/M 

=M 
DO 140 I=2,10-1 
assume that if T(I1). less 1.E-12 ,then no chem reac. 
TECUCL LIA = 2), GOrPO 140 
TRCEUIAG fh. 0799) GO TO 112 
GO/TO0-113 

CONTINUE 

CF1=Q*YF(1) 

CO1=FS*Q*YO(T) 

CFO=CF1+CO1 

CLA=0 .5*(CFO+ABS (CF1-CO1) ) 
CSM=CFO-CLA 

CF=CF1+T(1) 

CO=CO1+T (TI) 


CALL AS 
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CO 450, 121 
113 CONTINUE 
CF=T(1I)+(Q*YF(I)) 
CO=T (I) +(FS*Q*YO(I)) 
TB=T(I) 
LEGD CL en roe) CaNTo a9 
CALL CHEM 
T(I)=TB+REAC 
LF((CF-T(1)).LT.ERRO.OR.(CO-T(I)).LT.ERRO) THEN 
T(I)=TB 
GO TO 112 
END IF 
GO MEO? 
115 CONTINUE 
DO 120 J=1,MA 
TA=T (I) 
CALL CHEM 
AR=REAC 
T(I)=TA+(0.5%*AR) 
CALL CHEM 
BR=REAC 
T(1)=TA+(0.5*BR) 
CALL CHEM 
CR=REAC 
T(1)=TA+CR 
CALL CHEM 


T(1)=TA+( (AR+2 .*BR+2 .*CR+REAC) /6. ) 
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TEGCCF-T(1)).LT.ERROL.OR. (CO-T(1)).LT.ERRO1) THEN 
T(1I)=TB 
GO TO 112 
END IF 
120 CONTINUE 
121 CONTINUE 
pr kele) iC Fe IGE) s) /0 
YO(I)=(CO-T(L) )/(FS*Q) 
140 CONTINUE 
solve the temperature and species equation 
CALL, SY 
LECLELAG, EQ...2.) .GO «TO 150 
GO TO 160 
150 CONTINUE 
IF(INOWR.EQ.1) GO TO 152 
WRITE(10,*)'IFLAG,N,DT=’ , IFLAG,N,DT 
152 CONTINUE 
IFLAG=0 
DT=DT/2. 
NDUP=NDUP*2 
NHS=NHS*2 
NTD1=N+20 
MA=0 
DOwi55: I=1 , 10 
T(1)=TOD(T) 
YF(1)=YFOD(T) 


YO(I)=YOD(I) 
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155 CONTINUE 
GO TO 110 
160 CONTINUE 
check if there is ignition 
IF(N.LT.NIG) GO TO 200 
GO TO 210 
200 N=N+1 
TIME=TIME+DT 
GO TO 110 
210 CONTINUE 
IF(NIGS.EQ.1) GO TO 217 
LE(N.LE.NIGT) GO TO 213 
LEG GOLGN YW. ET20 45) GO TOR 2) 1 
IFGCTIME-TIMETD) .CT.TRES..AND.T(LIGN) .LT.0,2)) GO, TOswae 
IF((N-NIGT) .GT.4500.AND.T(IIGN) .GT.0.4) GO TO 400 
GO"TO.<213 
not dump data, that pick the data at N=N+l 
211 CONTINUE 
NHS=NHSTD 
NIGT=N+NHS 
N=N+1 
DT=DTTD 
NDUP=NDUPTD 
TIME=TIME+DT 
INOWR=0 
NCH=NCH/40 


GO TO 110 
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213 CONTINUE 
IF((N-220) .GT.NID)DT=DT3 
DO 214 I=2,1G 
Feat iy. Givi 0) THEN 
NIGS=1 
II=I 
WRITE GLOL*) pilot ign™at n,t,., l= ,N,TIME,1,1(f) 
TIMETDS=TIMETD-DT1 
WRITE(10,*) ‘NID, TIMETD=’ ,NTID, TIMETDS 
Micrel), Yr esOrarat i= TIDC LL) YEID( IL); YOID(IL) 
Peer iMETDS ; 11 
WhrieGini-) I,YF.YO,at ts lG=" -TIDCIG) ,YFID(1IG) , YOTDCIG) 
1+, TIMETDS 1G 
WRITE (1S5oF plot ign at n,t,i,I=’ ,N,TIME,1,T(I) 
WRITE(13,*) ‘NTD, TIMETD=' ,NTD, TIMETDS 
Hei tee Pareir yO ,at t,l=sTID(il), YFIDCIT), YOTD(IL) 
Peeer imei s , UL 
Weir eClo, se th Oy ateunhG= ,11D( 1G) YFIDCIG) , YOTDCLG) 
IPL IMETDS ,.1G 
END IF 
214 CONTINUE 
DO 216 I=IG+IGN+1,IN 
TPCT CL): GT 709.99)" THEN 
NIGS=1 
II=I 
WRITE(10,*)‘pilot ign at n,t,i,T=',N,TIME,1,T(I) 


TIMETDS=TIMETD-DT1 


216 


217 


220 
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WRITE(10,*)‘NTD, TIMETD=' ,NTD, TIMETDS 
WRITE(LO,*)'T,YF,YO,at t,I=* ,TID( LL) axhipe ll) | yOoup Gin 
, TIMETDS , II 
WRITE(10,%)’T,YF,¥O, at) t,1IG=' ,TTD(IG), .YFID( IG). YOTD< Te) 
PLIMETDSEUG 
WRETE(13).*)piloteten atun, t, 1, =A NeliMbe eae 
WRITE(13,*) ‘NTD, TIMETD=' ,NTD, TIMETDS 
WRITE(135*)‘T,YF,YO,at t,I=" ,TID(IL) , YeIDGLD) .yORDGL 
) PIMETDS,, TL 
WRITE(13.,.*)/T, YF,XO, at.t,1G=" , TTD( IG), YEED GIG), YOT Deis 
, LIMETDS, 1G 
END _ IF 
CONTINUE 
CONTINUE 
IF((N-400) .GT.NTD) DT=DT3 
NB=(N/NDUP) *NDUP 
IF(NB.EQ.N) CALL DUMP 
IF(NB.EQ.N) WRITE(10,*)‘N,DT=' ,N,DT | 
PECNIGSREQ > 1) .GORTO. 232 
turn on the ignition source 
IF(N.EQ.NIGT)GO TO 220 
COsTON23 2 
CONTINUE 
DO 225 I=1,10 
TTD(1I)=T(1) 
YFTD(1)=YF(1) 


YOTD(1I)=YO(I) 
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225 CONTINUE 
NTD=N+1 
NCH=NCH*40 
TIMETD=DT+TIME 
TeerDUFSgt 0.28 not turn on, Lene@sotiree 
DT=DT2 
DO 230 I=1,IGN 
T(IG+I)=1. 
230 CONTINUE 
232 CONTINUE 
N=N+1 
TIME=TIME+DT 
IF(N.GT.NMAX) GO TO 400 
IF(N.EQ.NTD1) THEN 
DT=DT2 
NDUP=NDUPTD 
NHS=NHSTD 
INOWR=1 
END IF 
GO TO 110 
400 CONTINUE 
REWIND (UNIT=12) 
WRITE(10,*)'N,t,DT=’ ,N, TIME, DT 
WRITE(12,*)'data for t,N=' ,TIMETD,NTD 
DO 410 I=2,10 
WRITE CL2, 1200) Te RCI) TID (1) ,YFID(1) , YOTDCT) 


410 CONTINUE 
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1100 FORMAT(’ 't, TFO2) 1G, IN" £8 .3) 1X ,9FS 42) 
1200 FORMAT ('}°D,X, TYR) YO= 4 5 I RUS soy.) 
NN=NMAX+450 
WRITE(12,*) ' NMAX+450=' , NN 
WRITE(13;*) ' NMAX+450=' ,NN 
IFCIFLAG.EQ.1) WRITE(10,*) ‘IFLAG,N=’ , IFLAG,N 
CLOSE(UNIT=10) 
CLOSE(UNIT=11) 
CLOSE(UNIT=12) 
CLOSE(UNIT=13) 
STOP 
END 
KAKKKKKA KAKA KEKE REE RE RERERERER ERE 
This subroutine is to calculate the solution of tri-diagonal matrix 


SUBROUTINE SY 


COMMON T(401) ,YF(401) ,YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, NDUM 
1 ,NIGS,DTOD,IM,IN,10,NODUF, TIME, DT, IFLAG, IDUMP,HL, NMAX, SMA 
LV be Dies DP2 DTS VCrYCOPAE  BT,.O; DAS FSUREAG. CLA CoMas le 
1 ,AFU(6) ,DF(6),BFC6), YFS 
DIMENSION DD(401) ,DD1(401) ,AA(401) , BB(401) 
DO 10 I=2,IM-1 
DD(1I)=D(2) 
AA(I)=A(2) 
BB(1I)=B(2) 
10 CONTINUE 


DD(IM)=D(3) 


20 


30 


40 


AA(IM)=A(3) 

BBC IM)=B(3) 

DO 20 I=IM+1,IN-1 
DD(1)=D(4) 

AA(1I)=A(4) 
BB(1)=B(4) 

CONTINUE 

DD(IN)=D(5) 

AA(IN)=A(5) 

BBC IN)=B(5) 

DO 30 I=IN+1,10-1 
DD(1)=D(6) 

AA(1)=A(6) 
BB(1)=B(6) 

CONTINUE 

solve the temp equ. 

DD(1)=1. 

AA(1)=0. 

DD(I0)=1. 

BB(10)=0. 

solve by Gauss elimination 

DD1(1)=DD(1) 

DO 40 I[=2,I10-1 
CC=BB(1)/DD1(I-1) 
DD1(1)=DD(I) - (CC*AA(I-1) ) 
T(1)=T(1) -(CC*T(I-1)) 


CONTINUE 
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50 


52 


54 
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DO 50 I=2,10-1 
TI=10-I+1 
T(II)=(T(II) -AA(II)*T(II+1) ) /DD1 (IT) 
IF(T(II).LT.1.E-22)T(II)=0. 
CONTINUE 
solve the fuel & Oxygen equ. 
DD(1)=DF(1) 
AA(1)=AFU(1) 
DD1(1)=DD(1) 
DO 52 I=2,IM-1 
DD(I)=DF(2) 
AA(I)=AFU(2) 
BB(I)=BF(2) 
CONTINUE 
DD(IM)=DF(3) 
AA(IM)=AFU(3) 
BB(IM)=BF(3) 
DO 54 I=IM+1,IN-1 
DD(I)=DF(4) 
AA(1)=AFU(4) 
BB(1)=BF(4) 
CONTINUE 
DD(IN)=DF(5) 
AA(IN)=AFU(5) 
BB(IN)=BF(5) 
DO 56 I=IN+1,I0-1 


DD(I)=DF(6) 


50 
AA(1)=AFU(6) 
BB(I)=BF(6) 
56 CONTINUE 
DO 60 I=2,10-1 
CC=BB(1)/DD1(I-1) 
DD1(1)=DD(I1) - (CC*AA(I-1) ) 
YF(1I)=YF(1) - (CC*YF(I-1)) 
YO(1I)=YO(T1) - (CC*YO(I-1)) 
60 CONTINUE 
DO 70 I=2,10 
II=IO-I+1 
YF(II)=(YF(I1) -AACII)*YF(II+1))/DD1 (IT) 
YO(II)=(YO(IIL) -AACII)*YO(II+1) )/DD1 (IT) 
DPOYE CID elie, B-22)YRC1I1)=0. 
IRCYOCLD) -LTI1LTE-22)YOC 11) =0. 
70 CONTINUE 
RETURN 
END 
KAKA AHA AKA AKAIKE RANI EEE 
This subroutine is to dump the data 
SUBROUTINE DUMP 
COMMON T(401) ,YF(401) ,YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, NDUM 
1 ,NIGS,DTOD,IM,IN,10,NODUF,TIME,DT,IFLAG, IDUMP,HL,NMAX, SMA 
1 ,1,DT1,DT2,DT3,CF,CO,AF,BT,Q,DA, FS,REAC, CLA, CSM, STP 
1 ,AFU(6) ,DF(6) ,BF(6) ,YFS 


IA=I0/IDUMP 


co. 
HX=0 . 8299 
NDUM=NDUM+1 
IF(NDUM.GT.1.5) GO TO 20 
TEQNIGS jEO M1) GOO 
NDUM=0 
REWIND (UNIT=12) 
5 CONTINUE 
IF(NIGS.EQ.1) WRITE(10,*)‘Pilot ignition at N,T=’ ,N, TIME 
IF(NIGS.EQ.1) WRITE(13,*)’Pilot ignition at N,T=’ ,NfTIME 
WRITE(12,*) ' N,TIME,DT=’ ,N,TIME,DT 
DO 10 I=2,10,IA 
HEX=HXACT CT) -T C1 - LO / OCI A) 
WRITE(12,2000)1,XX(1I) ,T(1) , YF(1I) , YO(1) , HFX 
10 CONTINUE 
GO TO 70 
20 CONTINUE 
IF(N.EQ.NMAX) GO TO 5 
check if there is diffusion flame 
IDUF=2 
TDUF=0. 
DOs 40) sl) sT0- 1. 
LFCTCD) .GT.TDUF) 240470430 
GO TO 40 
30 CONTINUE 
IDUF=I 
TDUF=T (IL) 


40 CONTINUE 


60 
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IFCTDUP oLT..0.28) THEN 
WRITE(10,*)‘'flame quenched at N,T=‘,N,TIME 
WRITE(13,*)'flame quenched at N,T=' ,N,TIME 
NDUM=0 
NIGS=0 
NODUF=0 
NMAX=N+25000 
REWIND (UNIT=12) 
CO4TO.5 
END (IE 
IF(((NDUM/20)*20) .EQ.NDUM) GO TO 60 
IF(NODUF.EQ.1) GO TO 70 
NODUF=0 
cc=0. 
IF(TDUF.GT.0.6.AND.YO(1).LT.0.001) THEN 
DT=DT3 
NMAX=N+4 500 
NODUF=1 
END IF 
CONTINUE 
IF(NODUF.EQ.0) GO TO 70 
WRETECLO's*) “These sis a diffusion flame at’ 
WRITE(10,*)‘1I,N,TIME=’ , IDUF,N, TIME 
WRITE(10,*)‘'SMA,STP=' ,SMA, STP 


WRITE(13,*)‘There is a diffusion flame at’ 


WRITEGL3 )*) °“1,N,TIME,DI=' , IDUF,N, TIME, DIE {SHR POI /GM2\S< 


I=IDUF 


Lo 


HFX=HX* (T(1) -T(1I-1))/(XX(1) -XX(I-1)) 
WRITE(1L3' 2000) XX D)., PCL)4 YRC Tay YO Gip, HES 
WRITE(12,*) ‘ N,TIME,DT=' ,N,TIME,DT 
DO 65 I=2,10,IA 
HEXMRHX* (TCL) SECT a1) )/ CAR (1) xX CI aD) 
WRITE (12 /2000)7 (| KX( I), TCL), YFCL)4 YOCL)S HEX 
65 CONTINUE 
70 CONTINUE 
2000 FORMAT.C’ &D, XT RYE SYO=9, 1S 1X (FLLI4a3 FIG, £9.28) 
RETURN 
END 
HAKAN KAKA KAKA KEKE KERR ERE KERR RE RRRREER 
This subroutine is to calculate the chemical reaction rate 
SUBROUTINE CHEM 
KKKAKKKKALAA KAKA KEKE KEK ERE KEKE KER ERE 
COMMON T(401) ,YF(401) ,YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, NDUM 
1 ,NIGS,DTOD,IM,IN,IO,NODUF, TIME,DT,IFLAG, IDUMP,HL,NMAX, SMA 
loi, DTL {DT2,DT3 {GF gCO SAFI BL,Q]DAsFS ,REAC,CLA,CSM, STP 
1 ,AFU(6) ,DF(6) ,BF(6) ,YFS 
Cl=1.-T(1) 
LECEGL CUT AUSE-20)° THEN 
GO TO 112 
END IF 
C2=BT*C1/(1.-AF*C1) 
IF (ABS (C2) .GT.30.) GO TO 112 
C3=1./EXP(ABS (G2) ) 


DPCCZVET.0.). C3=1,/C3 


154 
GO TO 114 
112 C3=0. 
PEC G2, CT: 07), LFLAC=1 
114 CONTINUE 
C4=C3*DT1/(FS*Q) 
REAC=DA* (CF-T(I) )*(CO-T(I) )*C4 


RETURN 


This subroutine is to calculate the chemical reaction rate near 
flame temperature 
SUBROUTINE AS 
KAKAKAKAKA KALLE KE KKKKEKEREREERERE 
COMMON T(401) ,YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, NDUM 
1 ,NIGS,DTOD,IM,IN,IO,NODUF, TIME, DT, IFLAG, IDUMP,HL, NMAX, SMA 
1 71, DELP) Crs CO" Ar, BT) DA, FS’, REAC, CLA, CSM, STP 
LD ., APUG.) DFCG), BF(6)% YFS 
Cl=1.-T(1) 
PATOL). LI 1 ke20 ) * THEN 
GO TO 112 
END IF 
C2=BT*C1/(1. -AF*C1) 
IF(ABS (C2) .GT.30.) GO TO 112 
C3=1./EXP(ABS (C2) ) 
TF (C2 -LT..0.) Co=175/C3 
GO TO 114 


thy bts a oe ee Ue 


55 


TF(GQeLTHOL). SFLAGST 


114 CONTINUE 


116 


120 


C4=DA*CLA*DT/ (FS*Q) 
LF(C4.LT:1.) GO TO 116 
IF(C3.GT.(50./C4)) THEN 
T(1)=T(1)+CSM 

co TO 120 

END IF 

CONTINUE 

C5=C4*C3 

C6=1./EXP(C5) 
T(1)=T(1)+CSM*(1.-C6) 
CONTINUE 

RETURN 

END 

this subroutine solves the solid equations numerically 


SUBROUTINE SFIN(SMA,STP,DT,TIME,HSF,CVF,TAMB,DNW,CNDA,APS,EPSL,N) 


COMMON/WDATA/WTD(451) ,WTDO(451) ,WAMD(451) , WAMDO(451) , WDN(451) 

1 ,WDNA(451) ,WDNC(451) , WHCP(451) ,WCND(451) ,WA(451) ,WD(451) , SIGM 
1 ,WD1(451) ,WDNO(451) ,WB(451) ,WR(451) , DNWF, DNF ,HCPA,HCPC,CNDC 

1 ,WDX,WTIME, EACT 

DOUBLE PRECISION CC1,CC2,CC3,WDT,WDX,WDN,WDNO 

DNW= Virgin wood density 

DNF= Final char density 

HCPA= Specific heat active wood 


HCPC= Specific heat of char 
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CNDA= Thermal condcutivity of active wood 
CNDC= Thermal conductivity of char 
EPSL= Emissivity 
HSF= Imposed heat flux 

CVF= h/EPSL ,h is convective heat transfer coefficient 
TAMB= Ambient temperature 

SIGM= 5.6696E-12 
APS= Ad, pre-exponential factor of wood 
NMAX= Maximum number of step 
WDX= Grid space of wood 

DT= Time step 

IO= number of grid space 

initial condition 

I0=350 

LEAN. Gia.) sGO.0b0.200 

change the unit of input data 

CNDA=0 .01*CNDA 

CNDC=0 .01*CNDC 

CVF=0 .O001*CVF/EPSL 
EACT=EACT*1000./1.987 

DO 100 I=1,10 

WID(1)=TAMB 

WIDO(1)=TAMB 

WDN(1)=DNW 

WDNO(1I)=DNW 

WAMD(1I)=0. 


WAMDO(I)=0. 


100 


160 


170 


200 


157 


CONTINUE 
SIGM=5 .6696E-12 

WDX2=WDX**2 

DNWF=DNW- DNF 

CONTINUE 

WTIME=TIME 

DO 170 I=1,10 

WDN(I)=WDNO(1) 

WTD(I)=WTDO(I) 

CONTINUE 

CONTINUE 

WDT=DT/0.55187 

IF((WIIME-TIME) .GT.1.E-12) GO TO 160 
WDNA(1)=DNW*(WDN(1) -DNF) /DNWF 

WDNC(1)=DNF* (DNW-WDN(1) ) /DNWF 
WHCP(1)=(WDNA(1)*HCPA+WDNC(1)*HCPC) /WDN(1) 
WCND(1)=(WDNA(1)*CNDA/DNW) +(WDNC(1)*CNDC/DNF) 
WA(1)=-WCND(1) /WDX2 

WB(1)=WA(1) 

WD(1)=(-WA(1) -WB(1))+(WDN(1)*WHCP(1) /WDT) 
WR(1)=WDN(1)*WHCP(1)*WTD(1) /WDT 

DO 210 I=2,I10 

WDNA(I)=DNW* (WDN(I) - DNF) /DNWF 

WDNG( I) =DNF* (DNW-WDN(Z) ) /DNWE 
WHCP(I)=(WDNA(1) *HCPA+WDNC (I) *HCPC) /WDN(1) 
WCND (1) =(WDNA (1) *CNDA/DNW)+(WDNC(I)*CNDC/DNF) 


WA(I)=-WCND(I) /WDX2 
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WB(I)=WA(I-1) 
WD(I)=(-WA(I) -WB(I) )+(WDN(I) *WHCP(I) /WDT) 
WR(I)=WDN(I)*WHCP(I)*WTD(I) /WDT 
210 CONTINUE 
WA(1)=1. 
WD(1)=-1. 
WB(1)=0. 
WA(I0)=0. 
WD(IO)=1. 
WB(I0)=-1. 
WR(1)=(-WDX*EPSL/WCND (1) )*(HSF-CVF*(WID(1) -TAMB) - SIGM*( (WID(1)**4) 
1 -(TAMB**4) )) 
where CVF=h/EPSL 
WR(I0)=0. 
solve the temperature equation 
DO 220 I=1,10 
WD1(I)=WD(T) 
220 CONTINUE 
DO 230 I=2,10 
CC=WB(I)/WD1(I-1) 
WD1(I)=WD1(I) - (CC*WA(I-1)) 
WR(I)=WR(I) - (CC*WR(I-1)) 
230 CONTINUE 
WR(10)=WR(1I0) /WD1(I0) 
WIDO(I0)=WTD(IO) 
WTD(I0)=WR(IO) 


DO 240 I=2,10 


240 


250 


260 
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II=10-I+1 
WR(II)=(WR(II) -WA(II)*WR(II+1) ) /WD1 (11) 
WTDO(II)=WTD(II) 
WID(II)=WR(II) 
CONTINUE 
WID(I0+1)=WTD(10) 
DO 250 I=1,10 
CC=EACT/WTD (I) 
CC1=APS /EXP(CC) 
APS=Ad 
CC2=CC1+(1./WDT) 
CC3=(WDN(I) /WDT)+CC1*DNF 
WDNO(I)=WDN(TI) 
WDN(I)=CC3/CC2 
CONTINUE 
DO 260 I=2,10 - 


II=I0-I+1 


WAMD (II )=WAMD(II+1) -WDX*(WDN( II) -WDNO(II))/WDT 


CONTINUE 
STP=(WTD(1) -TAMB) /(2230. -TAMB) 
SMA=WAMD (1) *3230. 

WTIME=TIME 

RETURN 


END 
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APPENDIX © 


PROGRAM LISTING FOR THE STEADY STATE DIFFUSION FLAME 
(FINITE DIFFERENCE METHOD AND ACTIVATION ENERGY ASYMPTOTICS) 


Finite difference numerical code 


PROGRAM FOR THE STEADY STATE DIFFUSION FLAME 
by: Lin-Shyang Tzeng 
solve T,YF,YO separately 
for steady state diffusing flame 
example for the input md of this program 
BEV0s7 eu e0Us7e) 0053/5072 051656 0.7 0 0.002 0°0 1.0 1.0 
OMS6089  6u1oe ze oro5b7 052510. 232 
POU me atte oioie soe el Oe Ol Los 10 
The result is writing on the file 'IGDUFLS.DAT' 
In the IGDUFLS.DAT, the last line SMA,SMAO are the fuel flow rate , 
SMA is the fuel flow rate where there is not a diffusion flame 
SMAO is the fuel flow rate where there is a diffusion flame 
PROGRAM ONEDSL 
COMMON T(401) , YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM,IN,10,NODUF,DX1,DX2 ,DX3, IFLAG, IDUMP 
2,PM1,NDST,YFS,RLE,AFU(6) ,BF(6) ,DF(6) ,ERR 
COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) , AA(401) 


1 ,BB(401) ,DDO(401) ,AAO(401) , BBO(401) 


160 


h616 \ 


DIMENSION TOD(401) , YFOD(401) , YOD(401) 
OPEN(UNIT=10, FILE=’ IGDUFLS .DAT’ , TYPE='NEW’ ) 

OPEN (UNIT=11, FILE='DAINS.DAT’ , TYPE='OLD’ ) 

OPEN (UNIT=12 , FILE='OUPTS ..DAT’' , TYPE=' NEW’ ) 

input data 
READ(11,*)DX1,DX2 ,DX3,SMA,STP, PARM, PM1, PSM, [PM,MX, YFS, RLE 
READ(11,*)AF,BT,Q,DA,FS,YOS 
READ(11,*)NMAX,NDST,IG,IM,IN,IO, IRA,RAD,NHS,NW, IDUMP,NCH 
IRA=1 including the radiation 

NDST=1 the density and the conductivity is not constant 
RLE ; Levise number 

PARM ; the ratio of the old and new data for the iteration 
PM1 ; maximum error for convergerne 

PSM ; error for the fuel flow rate of two 

NMAX ; maximum number for iteration 

LG 2einictral zt lame 7iocation 

IDUMP ; number of dumped data (IA=I0/IDUMP) 

RAD Planck mean absorption coeff Kp (/m) 

NCH,NHS undefined 

HL=(IM-1)*DX1+(IN- IM) *DX2+(1I0- IN) *DX3 

HL2=HL**2 

XF=1.-(LOG(0.058/YFS+1.)) /SMA 

IG=10*XF 

initial condition 

SMAO=SMA 

SMB=0. 


ERO=1.E-3 


102 


104 


106 


107 


sey 


coefficient of matrix and boundary cond. 
B(1)=0. 
A(1)=(-1.)/DX1 


AFU(1)=A(1)/RLE 


IFLAG=0 

XX(1)=0. 

DO 102 I=2,IM 
XX(1)=XX(1I-1)+DX1 

CONTINUE 

DO 104 I=IM+1,IN 

XX (1)=XX(1I-1)+DX2 

CONTINUE 

DO 106 I=IN+1,I0 
XX(1)=XX(1-1)+DxX3 
CONTINUE 

initial condition 
DO 107 I=1,1G-3 
T(1)=XX(1) /XX(IG) 
YF(1)=SMA*(1.-T(1)) 
YO(1I)=0. 
T(1)=T(1)+STP 
CONTINUE 

DO 108 I=IG+3,I0 
T(1)=(XX (10) -XX(T) ) /(XX(10) -XX(1G) ) 


YF(1I)=0. 
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YO(I)=YOS*(1.-T(I)) 
108 CONTINUE 
DO 109 I=1,7 
II=I-3 
TC(IG+II)=1. 
YF(IG+I1I)=YF(IG-4) 
YO(IG+II)=YO(1G+4) 
109 CONTINUE 
start the new iteration 
110 CONTINUE 
IF(SMA.LT.1.E-9) GO TO 400 
DO 111 I=1,10 
TOD(1I)=T(1) 
YFOD(1)=YF(T) 
YOD(1)=YO(1) 
111 CONTINUE 
LFCIPM: EQu1)eGOLTOSGO0 
ERR=0. 
calculate the old chemical reaction 
CALL CHEM 
DO 120 I=2,I10-1 
RTP=0. 
IF (IRA. EQ.1) THEN 
TP=T¢ 1 )*193204298 
DNST=0 .3528/TP 
RTP=RAD* (TP*¥*4)*1.576E-16/DNST 


END IF 
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RR(1I)=Q*YF(1I)*YO(1)*RR(I) -RTP 


Pee RRL) OUT .0.y-RRCL 0. 
120 CONTINUE 
CALL COEF 
IFC(IFLAG.EQ.1) GO TO 400 
solve the diffusion 
CALL SY 
DO 130 I=2,I10-1 
ER=ABS (T(1) -RR(I) ) 
IF(ER.GT.ERR) ERR=ER 
T(1I)=(1. -PARM) *T(1I)+PARM*RR(T) 
130 CONTINUE 
CALL CHEM 
CALL COEF 
IF(((N/2)*2).EQ.N) GO TO 142 
CALL SYF 
DO 135 I=1,I10-1 
ER=ABS (YF (1) -RR(I) ) 
IF(ER.GT.ERR) ERR=ER 
YF(I)=YF(1)*(1. -PARM)+PARM*RR(TI) 
135 CONTINUE 
CALL CHEM 
CALL COEF 
CALL SYO 
DO 140 I=1,I10-1 
ER=ABS (YO(I) -RR(T)) 


IF(ER.GT.ERR) ERR=ER 


140 


142 


145 


147 


148 
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YO(I)=YO(1I)*(1. -PARM)+PARM*RR (I) 
CONTINUE 

GO TO 148 

CONTINUE 

CALL SYO 

DO 145 I=1,10-1 

ER=ABS (YO(1) -RR(I)) 
IF(ER.GT.ERR) ERR=ER 
YO(I)=YO(1)*(1. -PARM)+PARM*RR(T) 

CONTINUE 

CALL CHEM 

CALL COEF 

CALL SYF 

DO 147 I=1,10-1 
ER=ABS (YF (I) -RR(I)) 
IF(ER.GT.ERR) ERR=ER 
YF(1I)=YF(1I)*(1. -PARM)+PARM*RR(1) 

CONTINUE 

CONTINUE 

IF(IFLAG.EQ.2) ERR=1. 

N=N+1 

IF(N.GT.NMAX) THEN 
WRITE(10,*) ‘not converge for N=’ ,N 
GO TO 400 

END IF 

IF(ERR.GT.ERO) GO TO 111 


ERO=1.E-4 
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DO 210 I=2,I0-1 
IF(T(I).GT.TMX) THEN 
TMX=T (I) 
LI=I 
END IF 
210 CONTINUE 
WRITE(10,*) ‘SMA, TMX, II, X=’, SMA, TMX, II, XX(II) 
IF(NW.EQ.0) GO TO 410 
IF(MX-1)220, 230,240 
220 CONTINUE 
IF(TMX.GT.0.27) THEN 
SMAO=SMA 
SMA=O . 5* (SMA+SMB) 
SMR=ABS (1. - (SMA/SMAO) ) 
IF(SMR.LT.PSM) THEN 
SMA=SMA*0. 5 
SMB=SMA 
END IF 
N=1 
GO TO 110 
END IF 
DO 225 I=1,10 
T(1)=TOD(I) 


YF(1)=YFOD(T) 


2a 


230 


tes 2 


234 


236 
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YO(1)=YOD(1) 
CONTINUE 
SMB=SMA 
SMA=0 . 5* (SMAO+SMA ) 
SMR=ABS (1. - (SMA/SMAO) ) 


IF(SMR.LT.PSM) GO TO 400 


‘Nel 


GO TO 111 
CONTINUE 
IF(TMX.GT.0.1) THEN 
DX1=DX1-0.0035 
DX2=DX2-0.00035 
DX3=DX3-0.0035 
A(1)=(-1.)/DX1 
DO 232 I=2,1M 
XX (I)=XX(I-1)+DX1 
CONTINUE 
DO 234 I=IM+1,IN 
XX(I)=XX(I-1)+DX2 
CONTINUE 
DO 236 I=IN+1,I0 
XX (1) =XX(1I-1)+DX3 
CONTINUE 
N=1 
GOMTORALO 
END IF 


GO TO 400 


240 


400 


410 
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CONTINUE 
WRLTECLO,*) 4ST (1) YRC@L)=4.T(1) , YF(1) 
IF(TMX.GT.0.27) THEN 
YFS=YFS-0.01 
N=1 
GO TO 110 
END IF 
CONTINUE 
SMA=SMA*1000 .0/3230.0 
SMAO=SMA0*1000 .0/3230.0 
IF(MX.EQ.0) WRITE(10,*)’minimum fuel flow rate (mg/cm2 sec) ,SMA 
1 ,SMAO=! , SMA, SMAO 
IF(MX.EQ.1) WRITE(10,*)’minimum quench dist.XX(I0) ,DX1,DX2 
Ly =i4 XX(IO}SDXIqDX2 
IF(MX.EQ.2) WRITE(10,*)‘min fuel con YFS=’ , YFS 
IF(IFLAG.EQ.2) GO TO 410 
IF(ERR.GT.0.01) THEN 
WRITE(10,*) ‘’ERR=’ , ERR 
GO TO 410 
END LE 
DO 410 I=2,10-1 
T(1)=TOD(1) 
YF(1)=YFOD(1) 
YO(1I)=YOD(T) 
CONTINUE 
CALL DUMP 


IF(IFLAG.EQ.2) WRITE(10,*)'IFLAG,N=’ , IFLAG,N 
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IF(IFLAG.EQ.1) WRITE(10,*)‘IFLAG,N=' , IFLAG,N 
CLOSE(UNIT=10) 
CLOSE(UNIT=11) 
CLOSE(UNIT=12) 


STOP 


This subroutine solve the tri-diagonal matrix 
SUBROUTINE SY 
KAKKAKAKA AKANE KERER EE KERERERRERERERE ER 
COMMON T(401) , YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM,IN,1I0,NODUF,DX1,DX2,DX3, IFLAG, IDUMP 
2. ,PM1,NDST;)YFS}RLE ,AFU(6) ,BFC6) ,DF(6)5 ERR 
COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) ,AA(401) 
1 ,BB(401) ,DDO(401) ,AAO(401) ,BBO(401) 
RR(1)=T(1) 
RR(10)=T(I0O) 
IF (NDSTIEQ? 1) “THEN 
DO 100 I=2,IM-1 
CUE=9007/ ChC1 )* 193 2e42965,) 
CDA=900./(T(1I+1)*1932.+298. ) 
CXT=1 . /(DX1*DX1*CDT) 
CXA=1./(DX1*DX1*CDA) 
SMX=SMA/DX1 
DD(1)=CXA+CXT -SMX 
AA(1)=SMX-CXA 


BB(1)=-CXT 
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100 CONTINUE 
CDT=900. /(T( IM) *1932.+298. ) 
CDA=900. /(T(IM+1)*1932.+298. ) 
CXT=1. / (DX2*DX1*CDT) 

CXA=1. /(DX2*DX2*CDA) 
SMX=SMA /DX2 
DD(IM)=CXA+CXT - SMX 
AA(IM)=SMX-CXA 
BB(IM)=-CXT 

DO 200 I=IM+1,IN-1 
CDT=900. /(T(I)*1932.+298. ) 
CDA=900. /(T(I1+1)*1932.+298. ) 
CXT=1 . /(DX2*DX2*CDT) 
CXA=1. /(DX2*DX2*CDA) 
SMX=SMA /DX2 
DD(1)=CXA+CXT - SMX 
AA(1)=SMX-CXA 

BB(1)=-CXT 

200 CONTINUE 
CDT=900. /(T(IN)*1932.+298.) 
CDA=900. /(T(IN+1)*1932.+298.) 
CXT=1. / (DX2*DX3*CDT) 

CXA=1. /(DX3*DX3*CDA) 
SMX=SMA/DX3 

DD (IN) =CXA+CXT - SMX 
AA(IN)=SMX-CXA 


BB(IN)=-CXT 


300 


10 


20 


DO 300 I=IN+1,I0-1 
CDT=900;/(T(1)*1932.+298.) 
CDA=900./(T(1I+1) *1932.+298.) 
CXT=1./(DX3*DX3*CDT) 


CXA=1./(DX3*DX3*CDA) 


SMX=SMA/DX3 


DD( 1) =CXA+CXT - SMX 


AA(I)=SMX-CXA 


BB(1)=-CXT 
CONTINUE 
GO TO 35 


ENDS LEG 


DO 10 I=2,IM-1 


DD(1)=D(2) 

AA(I)=A(2) 

BB(1I)=B(2) 
CONTINUE 
DD(IM)=D(3) 
AA(IM)=A(3) 


BB(IM)=B(3) 


DO 20 I=IM+1,IN-1 


DD(1)=D(4) 

AA(I)=A(4) 

BB(1)=B(4) 
CONTINUE 
DD(IN)=D(5) 


AACIN)=A(5) 


LAL 


BB(IN)=B(5) 

DO 30 I=IN+1,I0-1 
DD(I)=D(6) 

AA(1)=A(6) 
BB(1)=B(6) 
30 CONTINUE 
35 CONTINUE 

solve the temp equ. 

DD(1)=1. 

AA(1)=0. 

DD(1I0)=1. 

BB(I0O)=0. 

solve by Gauss elimination 

DD1(1)=DD(1) 

DO 40 I=2,10-1 
CC=BB(1)/DD1(I-1) 
DD1(1)=DD(1) - (CC*AA(I-1) ) 
RR(1I)=RR(1) - (CC*RR(I-1) ) 

40 CONTINUE 

DO 50 I=2,I10-1 
II=I0-I+1 
RR(II)=(RR(II) -AA(II)*RR(II+1) )/DD1 (IT) 
DEC(RRCEL) LIS Ce s12)) THEN 

RR(II)=PM1*RR(IL) 
IFLAG=2 
END IF 


50 CONTINUE 


ifs 
GO TO 80 
CONTINUE 
RETURN 


END 


KKK HAHA HAA KAKA HARARE KEKE KEKE REE 


This subroutine dump the data to a temporary file 


SUBROUTINE DUMP 


COMMON T(401) , YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM, IN, 10,NODUF,DX1,DX2 ,DX3, IFLAG, IDUMP 
2 ,PM1,NDST,YFS ,RLE,AFU(6) ,BF(6) ,DF(6) , ERR 

COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) ,AA(401) 
1 ,BB(401) ,DDO(401) ,AA0(401) , BBO(401) 

IA=10/IDUMP 

DO 10 I=1,10,IA 

WREPEC@LZ $2000) 17 XXCl) (LCL) VEG OCL) 

CONTINUE 

check if there is diffusion flame 

NODUF=0 

IDUF=2 

TDUF=0. 

DO 40 I=2,I0-1 

IF(T(I).GT.TDUF) GO TO 30 
GO TO 40 

CONTINUE 

IDUF=I 


TDUF=T (I) 
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40 CONTINUE 
WRITE(10,*)'max temp N,T,I=' ,N,TDUF, IDUF 
2000 FORMAT SL TYE .Y¥O«" (1S SURI SPT 4) 


RETURN 


shri kkkkkkob kkk bab kk kkk bak aka ak at abak ak ak abst obsbk 
This subroutine calculate the chemical reaction rate 


SUBROUTINE CHEM 


KKH KA KAKA KAKA KAKA KAKA IKKE KKK ERE ER REI 


COMMON T(401) , YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6),N, SMA 
1 ,RR(401) ,BT,AF,DA,FS, IM, IN, 10,NODUF, DX1,DX2,DX3, IFLAG, IDUMP 
2 ,PM1,NDST,YFS,RLE,AFU(6) ,BF(6) ,DF(6),ERR 
COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) , AA(401) 
1 ,BB(401) ,DDO(401) ,AAO(401) , BBO(401) 
DO 120 I=2,I10-1 
CleL3-TUD) 
C2=BT*C1/(1.-AF*C1) 
FPCABS(C2-CT00 ) Co MOM 
C3=1. /EXP(ABS (C2) ) 
EE (C2LT..0.0) (C3= cs 
GO TO 114 
112 c3=0. 
IF(C2.LT.0.) IFLAG=1 
IF(IFLAG.EQ.1) PRINT *,’I,N,T=',1,N,T(I) 
114 CONTINUE 
IF(NDST.EQ.1) THEN 


Density is not constant 


75 


DA1=DA*943 .4/(T(1)*1932.+298. ) 
GO “FOL ES 
END IF 
‘ DA1=DA 
115 CONTINUE 
RR(1)=C3*DAL 
120 CONTINUE 
RETURN 


END 


This subroutine calculate the coefficient of the tri-diagonal matrix 
SUBROUTINE COEF 

KKKKKKA KAKA KAA KA KAKA K AKAIKE HA EE 

COMMON T(401) ,YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM,IN,I0,NODUF, DX1,DX2, DX3 , [FLAG, IDUMP 
2 ,PM1,NDST,YFS,RLE,AFU(6) ,BF(6) ,DF(6) ,ERR 

COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) ,AA(401) 

1 ,BB(401) ,DDO(401) ,AAO(401) , BBO(401) 

D(1)=SMA+(1./DX1) 

A(2)=1.0*( (SMA/DX1) - (1. /(DX1**2) ) ) 

D(2)=0. - (SMA*1.0/DX1)+(2.*1.0/(DX1**2) ) 
B(2)=(-1.)*1.0/(DX1**2) 

A(3)=SMA*1 .0/DX2-2.*1.0/( (DX1+DX2)*DX2) 
D(3)=0.+2.*1.0/( (DX1+DX2) *DX2)+2 .*1.0/( (DX1+DX2)*DX1) 

1 -SMA*1.0/DX2 

B(3)=(-2.)*1.0/( (DX1+DX2) *DX1) 


A(4)=(SMA*1 .0/DX2) -1.0/(DX2**2) 
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D(4)=0. - (SMA*1.0/DX2)+(2.*1.0/(DX2**2) ) 

— B(4)=(-1.)*1.0/(DX2**2) 
A(5)=SMA*1.0/DX3-2.*1.0/( (DX2+DX3) *DX3) 
D(5)=0.+2.*1.0/( (DX2+DX3) *DX3)+2.*1.0/( (DX2+DX3) *DX2) 

1 -SMA*1.0/DX3 
B(5)=(-2.)*1.0/( (DX2+DX3) *DX2) 
A(6)=(SMA*1.0/DX3) -1.0/(DX3**2) 

D(6)=0. - (SMA*1.0/DX3)+(2.*1.0/(DX3**2) ) 
B(6)=(-1.)*1.0/(DX3**2) 
DF(1)=SMA+(1./(DX1*RLE) ) 
AFU(2)=(SMA/DX1) - (1. /( (DX1**2) *RLE) ) 
DF(2)=(2./( (DX1**2)*RLE) ) - (SMA/DX1) 
BF(2)=(-1.)/((DX1**2) *RLE) 
AFU(3)=SMA/DX2-2./( (DX1+DX2) *DX2*RLE) 
DF(3)=2./( (DX1+DX2) *DX2*RLE) +2. /( (DX1+DX2)*DX1 

1 *RLE) -SMA/DX2 
BF(3)=(-2.)/((DX1+DX2)*DX1*RLE) 
AFU(4)=(SMA/DX2) -1./((DX2**2) *RLE) 
DF(4)=2./((DX2**2)*RLE) - (SMA/DX2) 
BF(4)=(-1.)/((DX2**2) *RLE) 
AFU(5)=SMA/DX3-2./( (DX2+DX3) *DX3*RLE) 
DF(5)=2./((DX2+DX3) *DX3*RLE) +2. /( (DX2+DX3) *DX2 
1 *RLE) -SMA/DX3 
BF(S)=(-2.)/((DX2+DX3) *DX2*RLE) 
AFU(6)=SMA/DX3-1./( (DX3**2) *RLE) 
DF(6)=2./((DX3**2) *RLE) -SMA/DX3 


BF(6)=(-1.)/( (DX3**2) *RLE) 


ct ef 


RETURN 
END 
KEKKKAKK KAA KEKE KAA KE KEKE KALE ERE EE 
This subroutine solves the tri-diagonal matrix for fuel 
SUBROUTINE SYF 
KKKKKKKAKKAAK KAKA LAL NAKA AA NER ERE 
COMMON T(401) ,YF(401) ,YO(401) ,XX(401) ,A(6) ,B(6) ,D(6) ,N, SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM,IN,1I0,NODUF,DX1,DX2,DX3 , IFLAG, IDUMP 
2°); PM1,NDST,YFS,RLE,AFU(6),BF(6) , DF (6), ERR 
COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) , AA(401) 
- 1 ,BB(401) ,DDO(401) ,AAO(401) , BBO(401) 
RR(1)=SMA*YFS 
RR(1I0)=YF(10) 
IF(NDST.EQ.1) THEN 
DO 100 I=2,IM-1 
CDT=900./(T(1I)*1932.+298. ) 
CDA=900./(T(1I+1)*1932.+298. ) 
CXT=1 . /(DX1*DX1*CDT*RLE) 
CXA=1. /(DX1*DX1*CDA*RLE) 
SMX=SMA/DX1 
DDF(1)=CXA+CXT - SMX 
AAF(1I)=SMX-CXA 
BBF(1)=-CXT 
100 CONTINUE 
CDT=900./(T( IM) *1932.+298. ) 
CDA=900. /( TC IM+1)*1932.+298.) 


CXT=1 . /(DX2*DX1*CDT*RLE) 


200 


178 


CXA=1. /(DX2*DX2*CDA*RLE) 
SMX=SMA/DX2 

DDF (IM)=CXA+CXT - SMX 
AAF(IM)=SMX-CXA 
BBF(IM)=-CXT 

DO 200 I=IM+1,IN-1 

CDT=900. /(T(1)*1932.+298.) 
CDA=900. /(T(I+1)*1932.+298. ) 
CXT=1. / (DX2*DX2*CDT*RLE) 
CXA=1. /(DX2*DX2*CDA*RLE) 
SMX=SMA/DX2 

DDF (1) =CXA+CXT - SMX 
AAF(1)=SMX-CXA 

BBF(1)=-CXT 

CONTINUE 
CDT=900. /(T( IN) *1932.+298.) 
CDA=900. /(T(IN+1)*1932.+298. ) 
CXT=1. / (DX2*DX3*CDT*RLE) 
CXA=1. / (DX3*DX3*CDA*RLE) 
SMX=SMA/DX3 
DDF(IN)=CXA+CXT - SMX 
AAF(IN)=SMX-CXA 
BBF(IN)=-CXT 

DO 300 I=IN+1,I0-1 

CDT=900. /(T(I)*1932.+298. ) 
CDA=900. /(T(I+1)*1932.+298.) 


CXT=1./(DX3*DX3*CDT*RLE) 


300 


310 


10 


CXA=1. /(DX3*DX3*CDA*RLE) 
SMX=SMA/DX3 
DDF(I)=CXA+CXT - SMX 
AAF(1)=SMX-CXA 
BBF(1)=-CXT 

CONTINUE 
DO 310 I=2,I0-1 
DDF(I)=RR(I)*YO(1)+DDF(1I) 


CONTINUE 


CDA=900 ./(T(2)*1932.+298. ) 


DDF(1)=(1./(RLE*CDA*DX1) )+SMA 


AAF(1)=-1./(RLE*CDA*DX1) 

GO TO 32 

END IF 

DO 10 I=2,IM-1 
DDF(1)=RR(I)*YO(1)+DF(2) 
AAF(1)=AFU(2) 
BBF(1)=BF(2) 

CONTINUE 
DDF(IM)=RR(1)*YO(I)+DF(3) 

AAF(IM)=AFU(3) 
BBF(IM)=BF(3) 

DO 20 I=IM+1,IN-1 
DDF(I)=RR(1)*YO(1)+DF(4) 
AAF(I)=AFU(4) 


BBF(1)=BF(4) 


20 CONTINUE 


hao 


30 


32 


22 


40 
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DDF(IN)=RR(1I)*YO(1)+DF(5) 

AAF(IN)=AFU(5) 

BBF(IN)=BF(5) 

DO 30 I=IN+1,I0-1 
DDF(1I)=RR(1)*YO(1)+DF(6) 
AAF(1)=AFU(6) 
BBF(1)=BF(6) 

CONTINUE 

solve the fuel equ. 

DDF(1)=DF(1) 

AAF(1)=AFU(1) 

CONTINUE 

DO 35 I=2,10-1 
RR(1I)=0. 

CONTINUE 

DDF(10)=1. 

BBF(1I0)=0. 

solve by Gauss elimination 

DD1(1)=DDF(1) 

DO 4O I=2,I0-1 
CC=BBF(I)/DD1(I-1) 
DD1(1)=DDF(I) - (CC¥AAF(I-1)) 
RR(1)=RR(1) - (CC*RR(I-1) ) 

CONTINUE 

DO 50 I=2,10 
II=I0-I+1l 


RR(II)=(RRC(II) -AAF(II)*RR(II+1) ) /DD1 (IIL) 
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LFCRR(ALY SLT @- 2 E12 )) THEN 
RR(II)=PM1*RR(IT) 
IFLAG=2 
END IF 
50 CONTINUE 


RETURN 


This subroutine solves the tri-diagonal matrix for oxygen 
SUBROUTINE SYO 
KKKKKEKKAKA KEELER AER ERE ERER ERE RRR 
COMMON T(401) ,YF(401) , YO(401) ,XX(401) ,A(6) ,B(6) ,D(6),N,SMA 
1 ,RR(401) ,BT,AF,DA,FS,IM,IN,I0,NODUF, DX1,DX2,DX3, IFLAG, IDUMP 
2m, DMIs NDST YF SC RLESAFUCO) BECO) DEG] ERR 
COMMON DDF(401) ,AAF(401) , BBF(401) ,DD(401) ,DD1(401) ,AA(401) 
1 ,BB(401) ,DDO(401) ,AAO(401) , BBO(401) 
RR(1)=0. 
RR(1I0)=YO(I0) 
IF(NDST.EQ.1) THEN 
DO 100 I=2,IM-1 
SDI=J00M Cid) *1932.42082) 
CDA=900./(T(1I+1)*1932.+298.) 
CXT=1. /(DX1*DX1*CDT*RLE) 
CXA=1 . /(DX1*DX1*CDA*RLE) 
SMX=SMA/DX1 
DDO (1) =CXA+CXT - SMX 


AAO(1)=SMX-CXA 


100 


200 


BBO(1I)=-CXT 
CONTINUE 

CDT=900. /(T(IM)*1932.+298.) 
CDA=900. /(T(IM+1)*1932.+298.) 
CXT=1. /(DX2*DX1*CDT*RLE) 
CXA=1. /(DX2*DX2*CDA*RLE) 
SMX=SMA/DX2 
DDO(IM)=CXA+CXT - SMX 
AAO(IM)=SMX-CXA 

BBO(IM)=-CXT 

DO 200 I=IM+1,IN-1 

CDT=900. /(T(I)*1932.+298.) 
CDA=900. /(T(I+1)*1932.+298.) 
CXT=1. / (DX2*DX2*CDT*RLE) 
CXA=1. /(DX2*DX2*CDA*RLE) 
SMX=SMA/DX2 
DDO(I)=CXA+CXT - SMX 
AAO(1I)=SMX-CXA 

BBO(1L)=-CXT 

CONTINUE 
CDT=900. /(T( IN) *1932.+298. ) 
CDA=900. /(T(IN+1)*1932.+298. ) 
CXT=1. /(DX2*DX3*CDT*RLE) 
CXA=1. /(DX3*DX3*CDA*RLE) 
SMX=SMA/DX3 
DDO(IN)=CXA+CXT - SMX 


AAO (IN) =SMX-CXA 
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300 


310 


10 
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BBO(IN)=-CXT 
DO 300 I=IN+1,I0-1 
CDT=900. /(T(1)*1932.+298. ) 
CDA=900. /(T(I+1)*1932.+298. ) 
CXT=1. /(DX3*DX3*CDT*RLE) 
CXA=1. /(DX3*DX3*CDA*RLE) 
SMX=SMA/DX3 
DDO(I)=CXA+CXT- SMX 
AAO(I)=SMX-CXA 
BBO(I)=-CXT 
CONTINUE © 
DO 310 I=2,I0-1 
DDO(I)=DDO(I)+(RR(I)*YF(I) /EFS) 
CONTINUE 
CDA=900. /(T(2)*1932.+298. ) 
DDO(1)=(1. /(RLE*CDA*DX1) )+SMA 
AAO(1)=-1./(RLE*CDA*DX1) 
GOaTOss2 
END IF 
DO 10 I=2,IM-l 
DDO(I)=(RR(I)*YF(1) /FS)+DF(2) 
AAO(1I)=AFU(2) 
BBO(1)=BF(2) 
CONTINUE 
DDO(IM)=(RR(I)*YF(1I) /FS)+DF(3) 
AAO(IM)=AFU(3) 


BBO(IM)=BF(3) 


20 


30 


32 


39 
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DO 20 I=IM+1,IN-1 
DDO(1I)=(RR(1)*YF(1)/FS)+DF(4) 
AAO (1)=AFU(4) 

BBO(1)=BF(4) 

CONTINUE 

DDOCIN)=(RR(1I)*YF(1)/FS)+DF(5) 

AAO (CIN) =AFU(5) 

BBO(IN)=BF(5) 

DO 30 I=IN+1,I0-1 
DDO(1I)=(RR(1)*YF(1)/FS)+DF(6) 
AAO (1) =AFU(6) 

BBO(1I)=BF(6) 

CONTINUE 

solve the Oxygen equ. 

DDO(1)=DF(1) 

AAO (1)=AFU(1) 

CONTINUE 

DOV LZ OF 1 
RR(1I)=0. 

CONTINUE 

DDO(1I0)=1. 

BBO(1I0)=0. 

solve by Gauss elimination 

DD1(1)=DD0(1) 

DO 40 I=2,I10-1 
CC=BBO(1)/DD1(I-1) 


DD1(1)=DDO(T) - (CC*AAO(I-1)) 
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RR(I)=RR(I) - (CC#RR(I-1)) 
40 CONTINUE 
DO 50 I=2,10 
TI=I0-I+1 
RR(II)=(RR(II) -AAO(II)*RR(II+1) )/DD1(11) 
IF(RR(II).LT:.(-1.E-12)) THEN 
RR(II)=PM1*RR(II) 
LFLAG=2 
END IF 
50 CONTINUE 
RETURN 


END 


BKK HHA AAR EAR EAI EEE HEA TE IEA TERETE EE ETEK TK ae ak ae ae te ak ae ae ak ah ak 


Program for the calculation of extintion Damkohler number 


(Activation energy asymptotics method) 


This program is to calculate the minimum fuel flow rate by using the 
extinction Damkohler number derived by Linan | 

example for the input data, assume the fuel, is metnane 

SMA=0 .0002, YOX=0.232, HCM=1.5, TL=640.0, YFS=0.25, 
ACTV=48 .0, AEXF=2 .424E15, 

PROGRAM EXTIN 


FS=0.25 
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FS is the stoichiometry ratio of the fuel 
N=0 
PRINT *, ‘input SMA(g/cm2 sec) ,HCM(cm) , YOX,TL(k) , YFS’ 
READ(5,*)SMA,HCM, YOX,TL, YFS 
SMA is the fuel flow rate at the surface 
HCM is the boundary layer thickness 
YOX is the Oxygen mass fraction 
TL is the surface temperature, (absolute temperature k) 
YFS is the fuel mass fraction beneath the solid surface 
SMD=0 .000001 
PRINT *,’input activation energy, pre-exponedtial factor’ 
PRINT *,’kcal/mole,cm3/g sec’ 
READ(5 ,*)ACTV, AEXF 
ACTV is the activation energy of the gas fuel 
AEXF is the frequency factor of the gas fuel 
CONTINUE 
N=N+1 
IF(N.GT.1000) THEN 
PRINT *,’ NOT CONVERGE AT N=‘ ,N 
PRINT *, ‘REDAE,REDA’ ,REDAE,REDA 
GO TO 220 
END IF 
TFCSMA, LT). 1 E-1 2) 8COs TO. 200 
SMAD=SMA*3230. 
YFSO=YFS - (YFS+FS*YOX) /EXP(SMAD) 
TPND=35155.*YFSO 


STPL=TL/TPND 
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STPU=298 . /TPND 
DL=3230.*SMA 
ARFA=FS*YOX/YFSO 
BETA=STPL-STPU 
TACE=ACTV/(0.001983*TPND) 
TDEF=STPL+(ARFA-BETA) /(1.+ARFA) 
EBS=(TDEF**2) /TACE 
GAMA=1. -2.*(ARFA- BETA) /(1.+ARFA) 
ZETO=1.-(1./EXP(DL) ) 
DAMKO=4 . *0 .000374*AEXF*YFSO0*5 .61E-8/(0.323* (SMA**2) ) 
CN1=4 . *DAMKO* (ZETO**2 ) * (EBS**3 ) 
CN2=((1. - (ZETO*ARFA/(1.+ARFA) ) )**2)*((1.+ARFA) **2) 
CN3=1./EXP(TACE/TDEF) 
REDA=CN1*CN3 /CN2 
CN4=1.-GAMA 
REDAE=2 . 7183*(CN4- (CN4**2)+0.26* (CN4**3)+0 .055%* (CN4*%*4 ) ) 
IF(REDA.GT.REDAE) THEN 
SMA=SMA- SMD 
N=N+1 
GO TO 100 
END IF 
CONTINUE 
IF(( (REDAE-REDA) /(REDAE+REDA) ).GT.0.01) THEN 
SMA=SMA+SMD 
SMD=SMD*0 .1 
GO TO 100 


END IF 
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220 CONTINUE 
PRINT *,'SMA,HCM, YOX,REDA,REDAE=’ ,SMA,HCM, YOX,REDA, REDAE 
TDEFD=TDEF*TPND 
IF(N.EQ.1) PRINT *,‘'SMA too small,re-test again’ 
PRINT *,’diffusion flame temp.=’ , TDEFD 
STOP 


END 


Appendix D 


ESTIMATION OF SOLID PHASE PARAMETERS FROM 


PILOTED IGNITION EXPERIMENTAL DATA 


D.1. Experimental Data 


A series of piloted ignition experiments have been conducted by 
Abu-Zaid (1988). A combustion wind tunnel was used to perform these 
experiments in a well-controlled atmosphere under external radiation. 
This tunnel consists of three main sections, the inlet section, the test 
section, and the exhaust section; it is schematically shown in Figure 
D-l1. The crossectional area for air flow is 15cm x 8.5 cm, and the wood 
sample (14.5 cm long, 7.5 cm width and 3.75 cm high) is placed 30cm 
downstream of the leading edge. These wood samples were carefully 


instrumented by surface, bottom and in-depth thermocouples. 


During the experiments, the heaters were turned on at the desired 
radiant heat flux, and were allowed to warm up for about 15-20 minutes. 
The wood sample was then placed on the weighing table with its top 
surface flat along the tunnel base. A small natural gas pilot flame was 
used downstream of the tunnel to initiate piloted ignition. The output 
of gas analyzers, thermocouples and the load cell was stored and 


processed by a microcomputer in real time. 
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D.2. Estimation of Convective Heat Transfer Coeffeicient. 


The evaluation of the convective heat transfer coefficient during 
piloted ert tion is complicated, because (i) the surface temperature 
of the sample is continuously changing, and (ii) the heat losses from 
the surface do not occur primarily by convection. In fact, a major 
portion of the incident heat flux is lost by re-radiation. Thus, 
the use of correlations for the convective heat transfer coefficient 
obtained from either constant surface temperature or constant heat flux 
is questionable. Furthermore, there is a reduction in convective heat 
transfer due to blowing effects from the evolution of pyrolysis gases, 


although these are expected to be negligible for piloted ignition. 


In view of the above complications only an approximate value of 
convective heat transfer can be obtained from the Nusselt number 
correlations available in the literature. Also it seems more appropriate 
to use the constant surface temperature condition because the surface 
temperature at piloted ignition is nearly constant (since and the rate 
of change of surface temperature is small when it is in the neighborhood 
of the piloted tiniciba temperature). The estimation of the convective 
heat transfer coefficient is now decor t hed below. The mass flow rate in 
the wind tunnel during the piloted ignition experiments was about 
0.4625g/sec. Assuming the density of air in the wind tunnel to be 
0.8585 Kem. gives a free steam flow velocity of 4.22 cm/sec. Other 
properties of air (at 650 °K) are taken as: v = 58.510". merece 
aN (thermal conductivity of air) = 0.0495 wa K and Pr = 0.7. Also, 


the Nusselt number expression for laminar flow is used, viz., 


WAS) Ne yee 
Nuer0a32 po pea 
be 1G xX 


Since the leading edge of the thermal boundary layer is about 

30 cm downstream of the (actual) physical leading edge, Re. at the 
begining of the wood sample is Re(x=30cm)= ULx/v = 337.7Hence, 

the convective heat transfer coefficient is h(x=30 cm) = 0.8927 

pee The leading edge of the thermal boundary layer at the end of 
wood is about 44.5 cm. Thus, Re(x=44.5 cm) = 503.8 , and h(x=44.5 cm) = 
OEY poze! W/m K. Hence, the average convective heat transfer coefficient 

is taken as = 0.812 byte This value is used to calculate the surface 


temperature from the analytical (integral) solution. 


Although this estimate of the convective heat transfer 
coefficient is approximate, the convective heat transfer itself is not 
very important because the major heat loss from the surface occurs by 
radiation. The effect of h on the prediction of thermal conductivity 


by inverse calculations will be discussed later. 


D.3. Comparison of the Analytical Solution and Experimental Data 


In order to predict the experimental heat and mass transfer, a 
least square curve fitting program has been written to determine the 


optimal thermal conductivity. 


bag 


D.3.1 Determination of Thermal Properties of Wood 


In view of equation (5-17), if the convective heat transfer 
coefficient and the external heat flux have been determined, then the 
quantity inside the brackets is fixed. (say, equal to BX), and the 
constant K can be calculated by a least-square fit of the experimental 


values of-t'andnBXx. 


By using the equation (5-25) of Beck and Arnold (1977), the 


optimal value of K is given by 
2 
K = ut, BX, /Z(BX; ) : 


A computer program, SQFT, see Appendix E was written to calculate the 
optimal thermal properties from the measured surface temperature. 
Figure D-2 shows a comparison of the experimental and the 
least-square-error fitted surface temperature. The continuous line 

is the experimental data, and the dashed line is the least square fit. 
The thermal conductivity of the wood calculated by this procedure is 
listed in Table D-1. The average value of the thermal conductivity is 
0.25 W/mK, which is aici the range of literature-listed values. (the 
thermal conductivities for various woods, determined by Atreya (1983), 


range between 0.135 and 0.263 W/mK. (Assume rx, is a constant). 


TEMPERATURE (deg C) 


TEMPERATURE (deg C) 
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FIGURE D-2 Comparison of measured and predicted surface 


temperature: 


135 


Table D-1l. Thermal conductivity calculated from measured surface 


temperature. 


heat flux W/cm thermal conductivity W/mK 


D.3.2 The Effect of h on the Prediction of Thermal Conductivity of 


Wood. 


Since the major heat loss from the surface at piloted ignition 
is due to radiation, small changes in the convective heat transfer 
coefficient may not have a significant effect on the determination of 
the thermal conductivity of wood. Table D-2 shows the comparison of 
the thermal conductivity of wood calculated by least-square surface 
temperature curve fitting for different convective heat transfer 
coefficients. It is seen that even a 100% increase or decrease in 
the heat transfer coefficient results in only a small error in the 


determination of the thermal conductivity (~2%). 
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Table D-2. Comparison of thermal conductivity for different h. 
cia 
2 
W/cm | 


Oa315 O73:15 0.309 
On299 0.301 DM295 


D.3.3 The Pre-Expermental Factor Calculation 


The pre-exponential factor for wood pyrolysis is also not well 


known. In view of the fuel production rate equation, let ye = 155 A =b, 


and 
* * * 
(1-6 ) eo . A EXP(-E (1-1/T, ) 
x. = ——T Bee leas) ies > 
* *2 
x SE - = zs 


By using equation (5-25) Beck and Arnold (1977), the approperiate pre- 


exponential factor can be estimated from the equation 
* 2 
A = 27 .X./2%.., 
ts, i/ i 


where the thermal conductivity was assumed equal to 0.245 W/mK, 


uge 
the activation energy equal to 30 kcal/gm-mole, and the char yield § 
equal to 0.25. A computer program, SQFM, was written to calculate the 
optimal pre-exponential factor and the theoretical fuel evolution rate 


by using equation (5-18). The numerical code is listed in Appendix E. 


Figure D-3 compares the experimental and theoretical fuel flow 
rates. The continuous lines are the experimentally measured data, while 
the dashed lines are the fuel evolution rate calculated from the 
analytical solution. (The analytical solution for the pyrolysis of wood 
shows the same trend as the numerical calculation of pyrolysis at the 
very beginning). The shapes of the fuel evolution rate curves for the 
experimental data are different than the analytical solution curves. The 
experimental fuel evolution rate rises very quickly, then falls a bit; 
the drop to zero of the experimental data might occur because of the 
measurement error. Most of the gas flow out of the wood in the early 
stages of pyrolysis is water vapor. Wood contains much water (even "dry" 
wood). The experimental dry wood is only arbitrarily defined (heated 
at 105°C for 8 hours and allowed to cool over anhydrous CaSO,). The 
sharp rise for the fuel evolution rate at the beginning may occur 


because of the evaporation of the moisture contained in the wood. 


The pre-exponential factors calculated by the computer program 


are listed in Table D-3. 
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Table D-3. Pre-exponential factor for the pyrolysis, obtained 


by the best fitting of the experimental data. 


2 
heat flux W/cm pre-exponential factor (sec-1) 


average A = 1.296E8 (sec-1) 


The average pre-exponential factor is within the range of literature 


7 8 
values 6.0*10 ~7.5*10 (data from Tinney (1965)). 
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APPENDIX E 


PROGRAM LIST FOR THE PARAMETER ESTIMATION OF WOOD 


Program for the surface temperature 


(least square fitting for optimum heat flux and thermal cond.) 


This program least square fits the surface temperature of wood 
to estimate the convective heat transfer coefficient and thermal 
conductivity of wood 

PROGRAM SQFT 

DIMENSION T1I(800) , TA(800) ,CSS(800) , SSQTP(80) 

OPEN (UNIT=11, FILE='SQIN.DAT’ , TYPE=’OLD’ ) 
OPEN(UNIT=12 , FILE='SQOU.DAT' , TYPE='NEW’ ) 

PRINT) **) inpue conv, |, DNST EPSL in W/cem2,W/m2k,k,¢/cm35’ 

READ (5 ,*) CVF, TAMB,DNS, EPSL 

PRINT *,‘Input higher and lower conv Rads heat flux’ 
READ(5 , *)HSFH,HSFL 

SIGM=5 .6696E-12 

RLNTH=0 .126*TAMB/(HSF*85. ) 

CVF=CVF/EPSL 

HH=HH/EPSL 

HL=HL/EPSL 

HSF=HSF-(0.05/EPSL) 


HSF=HSFH 


200 


20 


201 

HSFD=(HSFH-HSFL) /80. 

PRINT *, ‘Input # of data’ 

READ(5,*)NDATA 

DO 20 I=1,NDATA 

READ(11,*)TI(I) , TACL) 

CONTINUE 
Least square fit of RKK then calculate sum of square of temp. 
CVF=HH 

CVFD=(HH-HL) /80. 
HSFM=0. 
RKKM=0. 

SSQTM=1 . E30 
HSFM is the minimum heat flux 

SSQTM is the minimum sum of square err SSQT 
RKKM is the minimum of RKK 
HSFM is dimensional, while RKKM is defined by Atreya’s thesis 
DO 400 Kl=1, 80 

CVF=CVF-CVFD 

HSF=HSF-HSFD 

AHTC=- (CVF*0 .0001+(4.*SIGM* (TAMB**3) ) ) 
SBC=-25 .*SIGM* (TAMB**2) /3. 
BTAB=(AHTC**2) - (4.*SBC*HSF) 

SSQ=0. 

SYX=0. 

DO 40 I=1,NDATA 

TEMPl=T-Tamb , Tamb=25. 


TEMP1=TA(I)-25. 
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AB1=HSF+AHTC*TEMP1+SBC* (TEMP1**2) 
CC=(TEMP1* (AHTC+2 . *SBC*TEMP1) /(BTAB*AB1) ) 
CCA=2 .*SBCTEMP1 
CCB=AHTC- (BTAB**0. 5) 
CCC=AHTC+(BTAB**0.5) 
CCD=(CCA+CCB) *CCC/( (CCA+CCC)*CCB) 
CC2=- (AHTC/(BTAB**1 .5) )*LOG (ABS (CCD) ) 
TIME1=(0 . 5*(TEMP1**2) /(AB1**2) -CC2-CC) 
SSQ=SSQ+TIME1**2 
SYX=SYX+TI (I) *TIME1 

40 CONTINUE 
RKK=SYX/SSQ 
TEMP=0. 
DTI=1. 
SSQT=0. 
DO 280 I=1,NDATA 
TIME=TI (I) 
TEMP1=TEMP+DTI 
DTI1=DTI 
IF(TIME.LE.1.E-7) GO TO 240 

100 CONTINUE 
AB1=HSF+AHTC*TEMP1+SBC* (TEMP1**2) 
CC=(TEMP1* (AHTC+2 .*SBC*TEMP1) /(BTAB*AB1) ) 
CCA=2 .*SBC*TEMP1 
CCB=AHTC- (BTAB**0 . 5) 
CCC=AHTC+(BTAB**0. 5) 


CCD=(CCA+CCB) *CCC/( (CCA+CCC) *CCB) 


203 
PRINT *, ‘CCD=' ,CCD 
CC2=- (AHTC/(BTAB**1.5) )*LOG(ABS (CCD) ) 
TIME1=RKK* (0. 5%* (TEMP1**2) /(AB1**2) -CC2-CC) 
ERR=(TIME1-TIME) /TIME 
IF(ABS(ERR).LE.1.E-3) GO TO 200 
IF(ERR.LT.O.) THEN 
TEMP=TEMP1 
TEMP1=TEMP+DTI1 
GO TO 100 
END IF 
DTI1=DTI1*0.5 
TEMP1=TEMP+DTI1 
GO TO 100 
200 CONTINUE 
SSQT=SSQT+(TEMP1-TA(I)+25.)**2 
280 CONTINUE 
HHF=CVF*EPSL 
IF(SSQT.LT.SSQTM) THEN 
SSQTM=SSQT 
HSFM=HSF 
RKKM=RKK 
END IF 
PRINT *, 'HSF,SSQT,RKK=' ,HSF, SSQT, RKK 
WRITE(12,1100)HSF, SSQT, RKK 
1100 FORMAT(1X, 'HSF,SSQT,RKK=' ,3F16.9) 
400 CONTINUE 


PRINT *, ‘Input HSF,RKK’ ,HSFM,RKKM 
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READ(5S,*)HSF,RKK 

CVF=CVF/EPSL 

AHTC=- (CVF*0.0001+(4.*SIGM* (TAMB**3) ) ) 

SBC=-25 .*SIGM*(TAMB**2) /3. 
BTAB=(AHTC**2) - (4.*SBC*HSF) 

TEMP=0. 

DTI=1. 

IF(CVF.LE.1.E-7) GO TO 480 

SSQT=0. 
REWIND (UNIT=12) 
DO 480 I=1,NDATA 

TIME=TI (1) 
TEMP1=TEMP+DTI 
DTI1=DTI 
420 CONTINUE 

AB1=HSF+AHTC*TEMP1+SBC* (TEMP1**2) 
CC=(TEMP1* (AHTC+2 .*SBC*TEMP1) /(BTAB*AB1) ) 
CCA=2 .*SBC*TEMP1 
CCB=AHTC - (BTAB**0 . 5) 
CCC=AHTC+(BTAB**0 . 5) 
CCD=(CCA+CCB) *CCC/((CCA+CCC)*CCB) 

CC2=- (AHTC/(BTAB**1.5) )*LOG(ABS (CCD) ) 
TIME1=RKK* (0. 5* (TEMP1**2) /(AB1**2) -CC2-CC) 
ERR=(TIME1-TIME) /TIME 

LF(ABS(ERR) .LE.1.E-3) GO TO 440 
IF(ERR.LT.O.) THEN 


TEMP=TEMP1 


440 


1200 


480 
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TEMP1=TEMP+DTI1 
GO TO 420 
END IF 
DTI1=DTI1*0.5 
TEMP1=TEMP+DTI1 
GO TO 420 
CONTINUE 
TEMP1=TEMP1+25. 
PRINT *,/TA(I),TEMP1=' , TA(I) , TEMP1 


WRITE(12,1200)TI(1I) , TEMP1 


FORMAT (2X, F10.5,1X,E12.5) 

CONTINUE 

STOP 

END 

KKKKA KAKA KARE KAKA KARE KKK KEKE KEKE REE EEE KEE RK KEK IK EK EEK EEE EKER 
Program listing of the least square fitting of pyrolysis of wood 
least square fitting for optimum pre-exponential factor 


PROGRAM SQFM 

DIMENSION TI(800) , FUMA(800) ,CSS(800) , SSQTP(80) 
OPEN(UNIT=11, FILE='SQIN.DAT’ , TYPE='OLD’ ) 
OPEN(UNIT=12 , FILE=’SQOU.DAT’ , TYPE=' NEW’ ) 


PRINT *,’input rads,T,DNST,EPSL,RTK,in W/cm2,k, g/cm3,W/mk’ 


206 
RTK is the thermal conductivity in W/mk 
READ(5,*)HSF,TAMB,DNS ,EPSL,RTK 
CVF=0.812 
RKK=0 . 66667*DNS*1.38*RTK*0.01/(EPSL**2) 
assume Cp = 0.33 cal/gk 
assume conv heat trans. coeff. = 0.812 W/cm2k 
SIGM=5 .6696E-12 
RLNTH=RTK*TAMB/(HSF*85. ) 
TD=(RLNTH**2 ) *DNS*1 .38/(RTK*0.01) 
CVF=CVF/EPSL 
HH=HH/EPSL 
HL=HL/EPSL 
PRINT *,'Input # of data’ 
READ(5,*)NDATA 
1100 FORMAT(2X,F10.5,14X,E12.5) 
DO 20 I=1,NDATA 
READ(11,1100)TI(1I) , FUMA(T) 
20 CONTINUE 
least square fit for optimum pre-exponential factor 
CVF=HH 
HSFM is the minimum heat flux 
SSQTM is the minimum sum of square err SSQT 
RKKM is the minimum of RKK 
HSFM is dimensional, while RKKM is defined by Atreya’s thesis 
AHTC=- (CVF*0 .0001+(4.*SIGM* (TAMB**3) ) ) 
SBC=-25.*SIGM* (TAMB*%*2) /3. 


BTAB=(AHTC**2) - (4. *SBC*HSF) 


420 
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SSQ=0. 

SYX=0. 

TEMP=0. 

DTI=1. 

IF(CVF.LE.1.E-7) GO TO 480 

SSQT=0. 

REWIND (UNIT=12) 

DO 480 I=1,NDATA 

TIME=TI(I) 
TEMP1=TEMP+DTI1 

DTI1=DTI 

CONTINUE 
AB1=HSF+AHTC*TEMP1+SBC* (TEMP1**2) 
CC=(TEMP1* (AHTC+2 .*SBC*TEMP1) /(BTAB*AB1) ) 
CCA=2..*SBC*TEMP1 

CCB=AHTC- (BTAB**0 . 5) 
CCC=AHTC+(BTAB**0 . 5) 
CCD=(CCA+CCB) *CCC/( (CCA+CCC) *CCB) 
CC2=- (AHTC/(BTAB**1 .5) )*LOG(ABS (CCD) ) 
TIME1=RKK* (0. 5* (TEMP1**2) /(AB1**2) -CC2-CC) 
ERR=(TIME1-TIME) /TIME 
IF(ABS(ERR).LE.1.E-3) GO TO 440 
IF(ERR.LT.O.) THEN 

TEMP=TEMP1 

TEMP1=TEMP+DTI1 

GO TO 420 


END IF 
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DTI1=DTI1*0.5 
TEMP1=TEMP+DTI1 
GO TO 420 
440 CONTINUE 
TEMP=TEMP1/TAMB 
TEMPS=TEMP+1. 
AB2=AB1/HSF 
RMS=(0. 75* (TEMPS**2) /(50.*AB2*EXP(50./TEMPS))*(1.-(1./ 
1 ((TEMPS**2)*EXP(50.*(1.-(1./TEMPS))))))) 
SMAND=FUMA (1) *0.001*TD/(DNS*RLNTH) 
SMAND is the non-dimensional mass flow rate 
SSQ=SSQ+(RMS**2) 
SYX=SYX+RMS*SMAND 
480 CONTINUE 
APND=SYX/SSQ 
APS=APND/TD 
APND is the non-dimensional pre-exponential factor 
APS is dimensional in sec-l 
TEMP=0. 
DTI=1. 
DO 280 I=1,NDATA 
TIME=TI (I) 
TEMP1=TEMP+DTI 
DTI1=DTI 
100 CONTINUE 
AB1=HS F+AHTC*TEMP1+SBC* (TEMP1**2) 


CC=(TEMP1* (AHTC+2 .*SBC*TEMP1) /(BTAB*AB1) ) 
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CCA=2 .*SBC*TEMP1 

CCB=AHTC- (BTAB**0.5) 

CCC=AHTC+(BTAB*%*0 . 5) 
CCD=(CCA+CCB) *CCC/( (CCA+CCC) *CCB) 

CC2=- (AHTC/(BTAB**1.5) )*LOG(ABS (CCD) ) 
TIME1=RKK* (0 . 5* (TEMP1**2) /(AB1**2) -CC2-CC) 
ERR=(TIME1-TIME) /TIME 
IF(ABS (ERR) .LE.1.E-3) GO TO 200 

EE GERRYLT 30. )) THEN 

TEMP=TEMP1 

TEMP1=TEMP+DTI1 

GO TO 100 

END IF 

DTI1=DTI1*0.5 

TEMP1=TEMP+DTI1 

GO TO 100 

200 CONTINUE 

TEMP=TEMP1/TAMB 

TEMPS=TEMP+1. 

AB2=AB1/HSF 

RMS=(APND*0 . 75* (TEMPS**2) /(50.*AB2*EXP(50./TEMPS) )*(1.-(1./ 
L ( (CTEMPS*¥2 )*EXP (5.00% C1y =iC1 2ATEMPS))))))) 
RMS FX=RMS*DNS*RLNTH*1000. /TD 

PRINT *,'RMS,FUMA(I)mg/cm2 sec’ ,RMSFX, FUMA(I) 
WRITE(12,*)TI(1I) , RMSFX 

280 CONTINUE 


PRINT *,‘APS (sec-1),APND,TD’ ,APS,APND,TD 
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CLOSE(UNIT=11) 


CLOSE(UNIT=12) 


STOP 


END 


HKKKKKKKKK KKK EKKKKEKERKEKEEEREREERE 


2h) 


LIST OF REFERENCES 


Abu-Zaid, M. [1988], "Effect of Water on Piloted Ignition of 
Cellulosic Materials", Ph. D. Thesis, Michigan State University, E. 
Lansing, MI. 


Aly, S. L. and Hermance, C. E. [1981], "A Two-dimensioal Theory of 
Laminar Flame Quenching", Combustion and Flame 40, 173-185. 


Amos, B., and Fernandez-Pello, A. C. [1988], "Model of the Ignition 
and Flame Development on a Vaporizing Combustible Surface in a 
Stagnation Point Flow: Ignition by Vapor Fuel Radiation Absorption", 
Combust. Sci. and Tech., 62,331-343. 


Anderson," DiMA. , Tannehill, J. C., and Pletcher.) R2 Hej teca 
"Combustional Fluid Mechanics and Heat Transfer", New York: 
Hemisphere Publishing Co. 


Atreya, A. [1983], Pyrolysis, Ignition and Fire Spread on Horizontal 
Surfaces of Wood", Ph. D. Thesis, Harvard University, Cambridge, MA. 


Atreya, A., Carpentier, C., and Harkleroad, M.[{1986],"“Effects of 
Sample Orientation on Piloted Ignition and Flame Spread,"First 


International Symposium on Fire Safety Science, V. 1, pp. 97-109. 


Atreya, A. and Wichman, 1.S8.[1987],"Heat and Mass Transfer during 
Piloted Ignition of Cellulosic Solids,"2nd ASME-JSME Thermal 


Engineering Joint Conference, V. 1, pp. 433-440. 


Atreya, A., and Wichman, I. S. [1989], "Heat and Mass Transfer 
During Piloted Ignition of Cellulosic Solids", Journal of Heat 
Traster.) 1) le pp are 


Bamford, C., Crank, J., and Malan, D., [1946], “The Combustion of 
Woodiy4Proce* Cambridge: Phil) Soct; 7 Voll7 42, 1p, 1667 


Baum, H. R., and Corley, D. M. [1981], "Finite Reaction Rate and Lewis 
Number Effects on Diffusion Flames in Channels", Presented at the 
Eastern State Section of the Combustion Institute, Pittsburgh, PA. 


Beck, J. V.,; and Arnold; K:.\J. [1977], “Parameter Estimation in 
Engineering and Science", New York, John Wilery & Sons. . 


Blackshear, P. L., and Kanury, A. M. [1970], "On the Combustion of 
Wood I: A Scale Effect in the Pyrolysis of Solids", Combust. Sci. 
and lech<4) 2, sop a. ae 


Buckmaster, J. D. [1975], "A New Large Damkohler Number Theory of Fuel 
Droplet Burning", Combustion and Flame, 24, pp.79-88. 


Coffee, T.P., Kotlar, A. J., and Miller, M.°S. [1982)) "other Overaiy 


2L2 


Reaction Concept in Premixed, Laminar, Steady-State Flames", 
Combustion and Flame, 54, 155. 


Gandhi, P. D. and, Kanury,,AwMe (1986)—. “Criteria for, Spontaneous 
Ignition of Radiantly Heated Organic Solids", Comb. Sci. Tech., 50, 
2554 


Glassman, I. [1977], "Combustion", New York, Academic press. 


Goos, A. (1952), "The Thermal Decomposition of Wood", Wood Chemistry, 
2, ACS Nomograph Series No. 97, 2nd Ed., Reinhold Pub. Co., New York. 


Hshieh, of. s., anderichards,, C.,N.. [1989], “Factors Influeticing 
Chemisorption and Ignition of Wood Chars", Combustion and Flame, 76, 


pp. o7-G7., 


Tshizuka,. Sv, andalsuji. HH. [1981], “An Experimental Study of the 
Effect of Inert Gases on Extinction of Laminar Diffusion Flames, 
"Eighteenth Symposium (International) on Combustion, The Combustion 
institute, Pittsburen,, PA,. pp.695-/03. 


Kanury, A. M., and Blackshear, P. L. [1970A], “Some Considerations 
Pertaining to the Problem of Wood-Burning", Combust. Sci. and 
Tech. #)l eppaso da). . 


Kanury, A. M., and Blackshear, P. L. [{1970B], “On the Combustion 
of Wood II: The Influence of Internal Convection on the Transient 
Pyrolysiceei. cellulose”, Combust. Sci. and Tech., 2,'pp.5-9. 


Kanury, A. M. [1972A],"Ignition of Cellulosic Solids - A Review,"Fire 
Research Abstracts and Reviews, V. 14, pp. 24-52. 


Kanury, A. M. [1972B], "Thermal Decomposition Kinetics of Wood 
Pyrolysis", Combustion and Flame, 18, pp./75-83. 


Kanury, A. M. [1974], “Ignition of Cellulosic Solids, Minimum Surface 
Temperature Criterion", Combust. Sci. and Tech:.; 9, pp.171-172. 


Kanury, A. M. [1975], “Introduction to Combustion Phenomena", New 
York, Gordon and Beach. 


Kaneda hee Pet lecw HB. and Chaiken,. Raw fF. [19/77]. “Mathematical 
Model of Wood Pyrolysis Including Internal Forced Convection", 
Combustion and Flame, 29, pp.311-324. 


Kashiwagi, T. [1974], "A Radiation Ignition Model of a Solid Fuel", 
Gomb. Sci! iTechateaisipp.225. 


Kashiwagi, T. [1981],"Radiative Ignition Mechanism of Solid Fuels, 
“Fire Safety Journal, V. 3, pp. 185-200. 


Kashiwagi, T. [1982], “Effect of Sample Orientation on Radiative 
Ignition", Combustion and Flame, 44,223-245. 


2s: 


Kashiwagiy 'T. ; Ohlemiller, Tra Jigeand Werner, K. 32987), sfbifece or 
External Radiation Flux and Ambient Oxygen Concentration on 
Nonflaming Gasification Rate and Evolved Products of White Pine", 
Combustion and Flame, 69, pp.331-345. 


Kent, J. H., and Wildiams, F. Ape(i9/4),  "Extinceionvot lamas 
Diffusion Flames for Liquid Fuels". Fifteenth Symposium 
(International) on Combustion, The Combustion Institute, Pittsburgh, 
Pennsylvania, pp.315-325. 


Kindelan, M. and William, F. A. [1975], "Radiant Ignition of a 
Comstible Solid with Gas-Phase Exothermicity", Acta Astronautica, 
2 ppeyass 


Kung, H. C. [1972], “A Mathematical Model of Wood Pyrolysis", 
Combustion and Flame, 18, pp.185-195. | 


Kung, H. C. [1974], “The burning of Vertical Wooden Slabs", Fifteenth 
Symposium (International) on Combustion, The Combustion Institute, 
Pittsburg. pp.243-253. 


Law, C. K. [1975], “Asymptotic theory for Ignition and Extinction in 
Droplet Burning". Combustion and Flame, 24, pp.89-98. 


Linan, A. [1974], "The Asymptotic Structure of Counterflow Diffusion 
Flames for Large Activation Energies", Acta Astronautica, l, pp. 
1007-1039. 


Martin, S. [1965], " Diffusion-Controlled Ignition of Cellulosic 
Materials by Intensive Radiant Energy", Tenth Symposium 
(International) on Combustion, pp.877. 


Nurbakhsh, S. [1989], "Thermal Decomposition of Charring Materials", 
Ph.D. Thesis, Michigan State University, E. Lansing, MI. 


Parker, W. J. [1988], "Prediction of Heat Release Rate of Wood” Ph.D. 
Thesis, The George Washington University. 


Puri, I. K., cand Seshadri, K. -(1986])),."“Extirnction‘tof Diffusion ykkemes 
Burning Diluted Methane and Diluted Propane in Diluted Air", 
Combustion and Flame, 65, pp.137-150. 


Roberts, A. F. [1970], "A Review of Kinetic Data for the Pyrolysis of 
Wood and Related Substances", Combustion and Flame, 14, pp.261-272. 


Sauer, F. [1956], "The Charring of Wood During Exposure to Thermal 
Radiation: Correlation Analysis for Semi-infinite solids", Intrim, 
T. AFSWP-868 U.S.D.A. Forest Service. 


Schwenker, R., and Beck, L (1963), “Study of the Pyrolytic 
Decomposition of Cellulose by Gas Chromatography", J. Polymer Sci. 
PartuGiml Nov2.4 ppeaade 


Sibulkin, M., Kulkarni, A. K., and Annamali, K. [1982], "Burning on 


214 


a Vertical Fuel Surface with Finite Chemical Reaction Rate", 
Combustion and Flame, 44, pp.187. 


Simms,-D. "Li" [1960], “"Denition” of Cellulosic Materials by Radiation", 
Combustion and Flame, 4, pp.293-300. 


Simms, D. L., [1961], “Experiments on the Ignition of Cellulosic 
Materials by Thermal Radiation", Combustion and Flame, 5, pp.369-375. 


Simms, D. L., [1962], “Damage to Cellulosic Solids by Thermal 
Radiation", Combustion and flame, 6, pp.303-318. 


Simms, D. L., [1963], "On the Pilot Ignition of Wood by Radiation" 
Combustion and Flame, 7, pp.253. 


SLs OL ence laws Mm i lL9os |). “The Ignition of Wet and Dry Wood by 
Radiation", Combustion and Flame, 11, pp.377-388. 


Spa ding bel OO mnt oorarc. Engineering 25,,..264; Proc. Roy. 
Soc, London 245A,, 352 (1958). 


Spa lerneee Ue bem oo seni tl Trans: Roy: Soc. London, A249, 1. 


izeng wel eoe ewe cnman st. Sand Baum, Hs Ri([:1987 )"AxCombined 
Analytical-Numerical Solution Procedure for Chemically-Reacting 
Flows", Central States Meeting of the Combustion Institute, pp.426- 
431. 


Tzeng, L. S., Atreya, A., and Wichman, I. S. [1990],"A One-Dimensional 
Model of Piloted Ignition", Combustion and Flame, 80, pp.94-107. 
1990. 


Wang, C. Y. [1972], "Lecture Notes, Perturbation Methods", Department 
of Mathematics, National Taiwan University. 


Meeson,~taRa,ewetker, ave RL Marid*Sliepcevichys Cr-Mg91(1971),¢The 
Piloted Ignition of Wood by Thermal Radiation", Combustion and 
Flame, 16, pp.303-310. 


Westprook, Co Key eand Dryer, F. Le [1981], “Simplified Reaction 
Mechanisms for the Oxidation of Hydrocarbon Fuel in Flames", 
Combustion Science and Technology, Vol.27, pp.31-43. 


Westbrook, C. K., and Dryer, F. L. [1984], “Chemical Kinetic 
Modeling of Hydrocarbon Combustion", Prog. Energy Combust. Sci., 
vee pp 1. -oO7,. 


Williams, F. A. [1982], “Urban and Wildland Fire Phenomenology", Prog. 
Energy Combust. Sci., 8, pp.317-354. 


Williams, F. A. [1985],"Combustion Theory,"The Benjamin/Cummings 
Publishing Co., Second Edition. 


Wichman, I. S., and Atreya, A. [1987], "A Simplified Model for the 


Zig 


Pyrolysis of Charring Materials", Combustion and Flame, 68, pp.231- 
247% 


Wichman, I. S., "On the Use of Operator-Splitting Methods for the 
Equations of Combustion", Combustion and Flame, to appear (1990). 


NIST-114A U.S. DEPARTMENT OF COMMERCE |1. PUBLICATION OR REPORT NUMBER 
(REV. 3-90) NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY NIST-GCR-91-595 
BIBLIOGRAPHIC DATA SHEET 


3. PUBLICATION DATE 
August 1991 


4. TITLE AND SUBTITLE 


Theoretical Investigation of Piloted Ignition of Wood 


5. 


AUTHOR(S) 


Lin-Shyang Tzeng and Arvind Atreya 


6. 


PERFORMING ORGANIZATION (IF JOINT OR OTHER THAN NIST, SEE INSTRUCTIONS) 
Michigan State University 

Department of Mechanical Engineering 

_East Lansing, MI 48824 


SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (STREET, CITY, STATE, ZIP) 
U.S. Department of Commerce 


National Institute of Standards and Technology 
Gaithersburg, MD 20899 


7. CONTRACT/GRANT NUMBER 


NIST Grant No. 60NANB8DO0861L 
8. TYPE OF REPORT AND PERIOD COVERED 


9. 


10. SUPPLEMENTARY NOTES 


11. ABSTRACT (A 200-WORD OR LESS FACTUAL SUMMARY OF MOST SIGNIFICANT INFORMATION. 
LITERATURE SURVEY, MENTION IT HERE.) 


IF DOCUMENT INCLUDES A SIGNIFICANT BIBLIOGRAPHY OR 


A theoretical model for piloted ignition of a flame in the gas phase above a vaporizing or 
pyrolyzing solid has been developed. Using this model it has been found that (i) The 
postulated simplified governing equations adequately explain the pre-ignition flashes that 
are often observed experimentally; (ii) A rational criterion for positioning the pilot flame 
exists; (iii) The heat losses to the surface play an important role, indicating that the fuel 
evolution rate by itself is insufficient for predicting the onset of piloted ignition. In 

this investigation, a numerical integration scheme is developed that accounts for the often 
vastly different rates between chemical reaction and convection or diffusion processes in the 
equations of combustion theory. This new numerical scheme is found to be very efficient for 
the piloted ignition problem, which involves both pre-mixed and diffusion flames. Finally, a 
numerical model for piloted ignition of wood which includes transient solid-phase decomposition 
has been developed. It has been found that the activation energy for the combustion of the 
evolved fuel is 49 Kcal/mole. 


12. KEY WORDS (6 TO 12 ENTRIES; ALPHABETICAL ORDER; CAPITALIZE ONLY PROPER NAMES; AND SEPARATE KEY WORDS BY SEMICOLONS) 


diffusion flames; equations; pilot ignition; premixed flames; wood 


13. AVAILABILITY 14, NUMBER OF PRINTED PAGES 


UNUMITED 
a 
3 


FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE TO NATIONAL TECHNICAL INFORMATION SERVICE (NTIS). 
X 


_ ELECTRONIC FORM 


15. PRICE 


ORDER FROM SUPERINTENDENT OF DOCUMENTS, U.S. GOVERNMENT PRINTING OFFICE, 
WASHINGTON, DC 20402. 


ORDER FROM NATIONAL TECHNICAL INFORMATION SERVICE 


S), SPRINGFIELD, VA 22161. 


a~ Le Se Ul le GTC cml tle mpibieaeeiahdeen 


~—ib~ Ee Plein There J rh "RSA oa veces Oy - 
ae YOOSOMMOET ‘OA bearyehh ne 49 ayurirewt J 


os es Sa ee ae nyis 21 On oem Ps ates) 


BE? “7aae ATAG SIHAAROO 


i] 
» } a. 2a ! , > : e 
i ical _——* ee ee ee a= pent ree 
7 wn aet ion. Alm] jee : “ppeae ts 
ke Herp4EAwa hata! rg 6 mois sal Inuovadg sty 
-~ —_— — _ a + a gee Seat: EE 
ave r2A balvwIA bos geet 
a UST SONT EN FED Tan NT REED AD TOs 
é 
& mitagiinee it 
; . Via T OS VR Bei PES 
; } > oe 
ye JED, 
~4 
= - - Oe atta 
; vow Thy ) hk ATOR 
{Len 
t Paty eee | basoliq ser lok 
i : sooleVvsh obad gad 
la: pabs | taupe goles 
T oe P ~ \ f - : r war! “ 
15 on P 29 ER 
Lisveaa 
§ i 435.39 
al IBG 
y a af 
i r ry bee « te 
J ’ rag 5 iat A 
: f H a } . i } 29 2B ‘ 
.sinm\la 
To lhe (@ @QnSw Ph4 grAmAcee Gud ;D kad PRO Tao © CRAs Ta Tere aay 
boow jeemslii bexiwweyq ~oottings sestq jenostaupe 7s 
} ’ eae 
> 
v4 a ig Oy 
—_ —— te raw tm ee RED mR : 
tar { Th + Yar . ¢, 
. os : P F : or y 
' ci A as 
i i ‘a 
~ Bt OO EONS MONT ODOC 1600597 GET A OT BAAAEEF THe OO lant 
— —_ — od ~~) 


uate ge Pe eras Tes OD A) STMUmOS 8a ica 


j j - >. : 4 a) . Beh csres 
acai = ee (ret vey Mery Om : no pe te ra 
7 =e . ‘pou eae, - ; 


- : © i] ay 
a - 


